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We consider the classical response in a chaotic system. In contrast to behavior in integrable 
or almost integrable systems, the nonlinear classical response in a chaotic system vanishes at long 
times. The response also reveals certain features of collective resonances which do not correspond 
to any periodic classical trajectories. The convergence of the response is shown to hold due to 
the exponential time dependence of the stability matrix. The growing exponentials corresponding 
to strong instability do not inhibit the convergence. We calculate both linear and second-order 
response in one of the simplest chaotic systems: free classical motion on a surface of constant 
negative curvature. We demonstrate the relevance of the model for applications to spectroscopic 
experiments. 
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I. INTRODUCTION 

Time-resolved femtosecond spectroscopy constitutes a 
powerful tool that probes electronic and vibrational co- 
herent dynamics of complex molecular systems in con- 
densed phase P, 0, S 13, IE Hi- Spectroscopic tech- 
niques allow for direct measurements of the time- and 
frequency-domain optical response functions that carry 
detailed information on the underlying dynamical phe- 
nomena. Nonlinear spectroscopy contains additional in- 
formation on noncquilibrium processes absent in linear 
measurements. 

In complex systems, such as polyatomic molecules, 
vibrational dynamics of many degrees of freedom can 
be adequately described in the framework of classical 
mechanics. In particular, at energies corresponding to 
room temperatures the complexity often originates from 
heavy strongly anharmonic vibrational modes that can 
be treated classically. 

A number of studies have been devoted to the non- 
linear response in stable (integrable) dynamical systems 
0) B B [iH ■ Classical nonlinear response functions 
have been shown to diverge with time in the general case 
although in certain cases the divergence can be re- 
moved by thermal averaging using the canonical distri- 
bution 0, [1] . The latter effect can be viewed as a result 
of destructive interference of different paths (dephasing), 
while the divergence is associated with the rephasing. 

The divergence would pose a major problem for the 
whole concept of perturbative response. On the other 
hand, quantum response functions do not exhibit any 
divergences. However, application of the fully quantum 
description is often unmanageable in the systems of in- 
terest. It has been demonstrated that systematic h ex- 
pansion cures the divergence when the limit of vanishing 
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h is taken after obtaining the asymptotics of long times 
d, • Alternatively the divergence can be removed by 
adding noise that models interaction with a bath (or a 
solvent) [13]. 

Previous analytical studies of the nonlinear response 
were mostly limited to fully solvable models, which im- 
plies integrable dynamics in the classical limit. Integrable 
dynamics corresponds to an idealized picture for a realis- 
tic physical situation, since even weak perturbations gen- 
erally destroy integrability. However, the problem is that 
weak deviations from the integrability do not eliminate 
the unphysical behavior of the response functions. This 
follows from the fact that the divergence in integrable 
systems originates from quasiperiodic motions on invari- 
ant tori ■ According to the Kolmogorov- Arnold-Moser 
(KAM) theory, most invariant tori are not destroye d by 
small perturbations that break down integrability [13| . 
Therefore, the long-time behavior of the response in al- 
most integrable systems is similar to its integrable coun- 
terpart. 

Stable dynamics is typical for a situation close to equi- 
librium. At larger energies a generic situation corre- 
sponds to dynamical behavior with chaotic features p^ . 
Moreover, unstable (hyperbolic) dynamics may be more 
common due to the stability of chaos with respect to 
perturbations. It has been argued based on results of 
numerical analysis [Tsj that chaotic dynamics appears to 
observe the convergence of the classical response func- 
tions. In spite of apparent importance and to the best of 
our knowledge, the problem of the nonlinear response in 
strongly chaotic systems has never been addressed using 
analytical methods. We note that analytical calculations 
of the response are rarely feasible in nonintegrable sys- 
tems, whereas numerical simulations of chaotic dynamics 
are complicated by the exponential divergence of stability 
matrices [T5| . 

A strongly chaotic (mixing) system is characterize d by 
a special spectrum of the Liouville operator [l^, [13, fig . 
The spectrum consists of complex Ruelle-Pollicott (RP) 



2 



resonances that determine the asymptotic oscillations 
and decay of the correlations. This yields the linear re- 
sponse function directly related to the two-point correla- 
tion functions, in agreement with the fluctuation dissipa- 
tion theorem (FDT). Nonlinear response, however, turns 
out to be more involved still with the noticeable effect of 
the resonances. 

In this work we show that (i) the classical linear and 
nonlinear response of a chaotic system exhibits decay and 
oscillations as a function of times between the driving 
pulses, and (ii) the Fourier transform of 2D second-order 
response function reveals broad and asymmetric peaks 
that can be viewed as signatures of chaos in underlying 
dynamics. 

We first establish a general qualitative picture of lin- 
ear and nonlinear response in a classical system with 
hyperbolic chaotic dynamics and demonstrate that the 
classical response functions exponentially vanish at large 
times. To exemplify our arguments we further perform 
detailed analytical calculations of linear and second-order 
response functions for a free particle moving along a com- 
pact surface of constant negative curvature. The model, 
that constitutes a well-known example of classical chaos 
and a prototype for quantum chaos, allows for an exact 
solution due to strong dynamical symmetry. To the best 
of our knowledge, we present the first explicit analytical 
calculation of nonlinear response in a chaotic system. 

The manuscript is organized as follows. We begin with 
reviewing the classical response theory using the Liouville 
representation of classical mechanics. In Section IIIII we 
present the qualitative picture for the classical nonlinear 
response functions and argue that they show exponen- 
tial decay at long times. We also argue that a generic 
finite motion of the molecular system with potential in- 
teractions is similar to the free motion along an effective 
compact configuration space of negative curvature. In 
Section [TVl we consider the chaotic model of free motion 
along a Riemann surface of constant negative curvature. 
We demonstrate how to make use of strong dynamical 
symmetry, and reduce the model to a much simpler prob- 
lem of ID dynamics on a circle. In Section |V] we present 
our calculation of the linear and second-order response 
functions. This calculation constitutes the main techni- 
cal result of the paper. The results of numerical calcula- 
tions of 2D spectroscopic signals are presented in Section 
IVIl All necessary technical details are summarized in the 
Appendices. 



II. CLASSICAL RESPONSE IN LIOUVILLE 
SPACE 

In this section we review a general formalism of the 
classical response theory, using the language of proba- 
bility distributions (e.g. see Ref. HI). We adopt the 
Liouville picture of classical mechanics where the system 
state that evolves in time is described by a distribution 
p(?7, i), whereas the observables are represented by func- 



tions in phase space. This language is advantageous in 
a chaotic case, since mixing (e.g., hyperbolic) dynamics 
is characterized by strong instability in trajectory space, 
and individual trajectories do not provide meaningful in- 
formation on the system relaxation. 

Consider classical Hamilton dynamics that occurs in 
phase space M equipped with a Poisson bracket. Evo- 
lution of an externally driven system is described by the 
classical Liouville equation with a time-dependent Hamil- 
tonian Hxit) 

'J^ = -Lrit)piv,t), (1) 
LT{t)... = {HT{t),...} , HT{t)=H-£{t)f , (2) 

where H{r]) is the time independent Hamiltonian of the 
undriven system, £{t) is the time dependent uniform 
driving field, and dipole moment (polarization) /(jy) de- 
scribes the coupling of the system to the driving field. 
This sign choice of in the Poisson bracket {/, g} cor- 
responds to a convention that the action of the Liou- 
ville operator Lt determines the phase space velocity as 
r) = LT{t)r}. 

In many cases /(j?) is also an observable, e.g. polariza- 
tion that creates spectroscopic signals, and the measured 
value is given by the time dependent average 

(/) = I dr,p{r,,t)f{rj). (3) 

The undriven system with Hamiltonian H is character- 
ized by equilibrium distribution po{ri) that depends only 
on the conserved energy H of the system, i.e. {H, po} = 
0. The average of / for the equilibrium distribution poiv) 
normally vanishes (in spectroscopic literature this is of- 
ten referred to as the absence of the permanent dipole): 

J dripo{v)f{v)^0- (4) 

The response theory is usually based on a perturbative 
expansion of the measured signal in powers of the driving 
field. The response functions can be naturally obtained 
via calculating the correction to the equilibrium density 
distribution due to the driving field. The full Liouville 
equation is repeatedly solved, this results in an expansion 
P — Po + Pi + ■ ■ ■ where pn represents the contribution of 
ri-th order in the field £. Therefore, we start with solving 
the equation 

{dt+L)p,^£it){f,po}, (5) 

to get the linear in £ correction to the density distribu- 
tion: 

t 

pi{t) = J dr^ £(ri)e-^(*--^){/, po} ■ (6) 



Here L is the Liouville operator 

Lp = {H,p} (7) 
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of the unperturbed dynamics that corresponds to the 
Hamiltonian H . The corrections to po in all orders of £ 
can be obtained by solving equations similar to Eq. ([5]). 
This leads to the natural representation of the observable 
quantity in the form of an expansion 

t 

(/) ^ j drip{v,t)f{rj) = I dnSi{t;n)£{n)+ (8) 



t T2 

J dT2 j dT, S2{t;T,,T2)£iT,)£iT2) + . . . , (9) 


which defines response functions S'„(t; . . . , ti). The re- 
sponse function of order n depends on n + 1 time vari- 
ables, it can be reduced to n time segments tn = t — 
r„, . . . ,ti = T2 — Ti in the case when the unperturbed 
dynamics has no explicit time dependence. For time 
segments representation, the linear and second-order re- 
sponse functions read: 

S^'Hh)^ j dvf{ri)e-^'Hf{v),Po}, (10) 

(11) 

These expressions for the first- and second-order response 
will be the starting point of our general qualitative anal- 
ysis of response functions in hyperbolic systems. They 
will be also utilized in our analytical calculations for a 
particular chaotic system where strong dynamical sym- 
metry of the phase space enables us to solve the problem 
completely. 

III. QUALITATIVE PICTURE OF RESPONSE 
IN A HYPERBOLIC DYNAMICAL SYSTEM 

A. Integrable and almost integrable dynamics 

Dynamics of a classical Hamiltonian system in the 
vicinity of an equilibrium (stationary) point can be ad- 
equately described by a system of uncoupled harmonic 
oscillators. This is achieved by expanding the classical 
Hamiltonian H up to second order in dynamical variables 
(the first-order terms vanish for an expansion around a 
stationary point). A harmonic system represents the sim- 
plest example of integrable dynamics. 

The evolution of integrable systems with n degrees 
of freedom is naturally described in terms of n canon- 
ical pairs of action and angle variables. Action vari- 
ables Cj with j = 1, . . . , 71 are first integrals of motion 
{H,Cj} = in involution, i.e. {ci,Cj} — 0, with the 
Hamiltonian H(r]) = H(ci, . . . , c„) depending on the in- 
tegrals of motion only. The corresponding phase space 
vector fields cj defined by cjf = {cj,f}, therefore, com- 
mute and describe motions along the tori with frequencies 



ujj = dH{ci, . . . ,Cn)/dcj that depend parametrically on 
the integrals of motion. Since for fixed values of the inte- 
grals of motion the angular velocities are constant, such 
dynamics is referred to as is conditionally periodic (or 
quasi-periodic) motion. For given values of action vari- 
ables determined by initial conditions, the trajectories 
lie on the corresponding n-dimensional torus which is a 
subspace of 2n-dimensional phase space. The motion on 
the torus is strictly periodic if the ratios of correspond- 
ing angular velocities are rational (resonant torus); if the 
angular velocities are incommensurate, the trajectories 
densely cover the torus (non-resonant torus). Harmonic 
dynamics represents the simplest particular case of inte- 
grable dynamics when the frequencies uji, i — I, . . . ,n are 
independent of the action variables Cj . 



The nonlinear response function for classical integrable 
dynamics have been shown to have power-like divergence 
at long times which can be eliminated by invoking a quan- 
tum description 0, d, H, [l^i [HI ■ Divergence of nonlinear 
response in integrable systems does not imply unphysical 
behavior by itself, since integrable dynamics is an ide- 
alization. In real systems quantum effects, interactions 
with the bath or irregular dynamics may provide a neces- 
sary regularization. While the first two phenomena have 
been discussed in the literature, we concentrate on the 
latter. A small, yet generic, perturbation of the system 
Hamiltonian destroys integrability. However such a per- 
turbation does not break down the linear divergence in 
the response functions. This follows from the qualitative 
picture of almost integrable dynamics established by the 
celebrated Kolmogorov-Arnold-Moser (KAM) perturba- 
tion theory [l3|. The KAM theory states that a suffi- 
ciently small perturbation does not destroy most nonres- 
onant tori, which means that in the invariant subspace 
of the entire phase space represented by the remaining 
distorted tori the dynamics preserves its quasiperiodic 
nature. Although motion inside the instability zones rep- 
resented by destroyed tori becomes chaotic, their relative 
measure in phase space is small for small perturbations 
to the integrable dynamics. In the simplest n = 2 case 
the zones of instability are confined between remaining 
invariant tori. 



It has been pointed out that the divergence or the 
classical response functions in integrable systems origi- 
nates from quasiperiodic nature of the underlying mo- 
tion. Combined with the picture of almost integrable 
dynamics established by KAM theory, this demonstrates 
the unphysical divergence of the response functions is 
stable with respect to at least weak deviations from in- 
tegrability. We will see that chaotic regions (instability 
zones) do not contribute to the divergence, therefore the 
diverging terms will decrease with increasing the devia- 
tion from the integrable situation due to decrease of the 
amount of surviving invariant tori. 
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B. Chaotic dynamics 

Integrable and almost integrable Hamiltonian dynam- 
ics considered in subsection IIII Al is typical for low ener- 
gies when the system is moving in the neighborhood of 
a stable stationary point (equilibrium) . Islands of stable 
dynamics completely vanish for larger deviations from 
the integrability caused by higher characteristic energies. 
The motion will be typically chaotic for higher temper- 
atures or in essentially nonequilibrium processes such as 
photoinduced molecular dynamics on the excited elec- 
tronic adiabatic surfaces. 

In this subsection we present a qualitative picture of 
nonlinear response in chaotic systems. Due to FDT, the 
stability matrices were found to affect the response func- 
tions starting with the second order [H, [l^. Therefore, 
the divergence of the nonlinear response functions in the 
integrable case can be attributed to the growth with time 
of certain stability matrix components. In the case of 
strong chaos there are components of the stability matri- 
ces that grow exponentially ~ e^' with time, A being the 
Lyapunov exponent. Numerical simulations have demon- 
strated that at least for some examples of chaotic dynam- 
ical systems the classical nonlinear response functions are 
free of unphysical divergences [31 ■ 

In this subsection we present qualitative arguments 
that rationalize physical exponentially decaying at long 
times behavior of classical nonlinear response functions 
in systems that exhibit strong enough chaotic behav- 
ior. More precisely we will consider mixing systems also 
known as A-flows, Smale or uniformly hyperbolic dynam- 
ical systems (for a nice overview without too many de- 
tails see Ref. |2CI) . In these systems at all points of the 
{2n — l)-dimensional energy shell, or all relevant points 
(the relevant points belong to the so-called nonwandcring 
subset) in the Smale's case, one can define unstable 
and n_ stable tangent directions with n-^- + n_ — 2n~2, 
so that a small deviation from a trajectory along the 
stable (unstable) directions decays exponentially in time 
for forward (backward) dynamics. The stable (unstable) 
directions can be locally integrated to obtain stable (un- 
stable) manifolds of dimension n_ («+). Note that n+ 
(n_) is the number of positive (negative) Lyapunov ex- 
ponents. Also note that in the Hamiltonian dynamics 
case the volume in phase space is conserved so that all 
Lyapunov exponents sum up to zero 

We will consider the simplest case of the lowest di- 
mension coordinate space dimension n = 2 that al- 
lows for chaotic dynamics. The isoenergetic shell is 3- 
dimensional, we have n+ — n_ = 1, and two nontrivial 
Lyapunov exponents ±A. 

A schematic picture of first- and second-order response 
formation in a hyperbolic chaotic system is shown in 
Fig. [T] Since chaotic dynamics described in terms of 




FIG. 1: Schematic picture of the cross-section of the phase 
space along the surface defined by the stable and unstable 
directions. Initial distribution of / is presented by two regions 
oppositely charged regions (dark red and light green) . As time 
elapses the distributions elongate along unstable directions 
and contract along stable ones. 



phase-space trajectories is extremely unstable, using the 
dual representation via evolution of distributions (Liou- 
ville representation) will be advantageous. We start with 
the linear response function S'-^^t). Since the equilib- 
rium distribution depends on the system energy only, we 
can recast the expression for the linear response [Eq. (jlOp ] 
in a form that actually represents the FDT 

S^^\t) = dtl^'\t) ^dtj drjfe-^*fdEPo ■ (13) 

Since J dxf — (no permanent dipole), we can consider 
the simplest set-up when the function / is represented 
by two small separate regions in the phase space where 
it adopts positive and negative values (an extension of 
our arguments to the general case is straightforward). 
By our assumption the linear sizes of the regions a <^ I 
are small compared to the linear size I of the compact 
phase space. Evolution during time t changes the shapes 
of the regions. In the case of hyperbolic dynamics for 
t 3> 1, the shape becomes similar to ribbon-like fettuc- 
cine: elongated along the unstable direction by a factor 
e^^, narrowed along the stable one by e~^* (we reiterate 
that A is the positive the Lyapunov exponent) and un- 
changed along the flow. The long ribbon becomes evenly 
folded in the entire accessible energy shell. The overlap 

of the distributions / and e"'"^ fdsPa in Eq. (fT3l) is repre- 
sented by a large number Ni(t) of disconnected regions. 
The dimensions of each region are typically a, ae~^* and 
a along the flow, in the stable and unstable direction, re- 
spectively. This gives the volume of the single 3D region 
in the energy shell vi{t) ~ a^e~^*. A typical distance 
d{t) between the fettuccine in the overlap regions can be 
easily estimated as d{t) ~ /■^a~^e~'^*. This estimates the 
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number of disconnected regions of the overlap as 

7Vi(t)^(//a)V*. (14) 

Since / assumes opposite signs in two initial regions, 
the distribution e~^^ fdEPo consists of two positively and 
negatively "charged" fettuccine. The cancellations result 
in a signal determined by a typical fluctuation propor- 
tional to -v/lVi, and the overlap integral attains a factor 
y/WiVi. This results in the decaying long-time asymp- 
totic of the linear response function S'(i)(i) - e'^^/"^. 

We are now in a position to consider the second-order 
response function S^^\ti,t2)- Propagating the observ- 
able / in Eq. (jlip backward in time we interpret the 
second-order response as an overlap 

s^^KhM) ^ dtj^^Hh,t.2) ^ dt, j dqu^^djf+ (15) 

of the distributions f-{ri) = exp{Lt2)f and ^^djf+{T]) 

with /+ = exp{—Lti)fdEPo and the vector field ^^dj = 
{/, •}. Similar to the case of linear response for t ^ 1, 
the shapes of both distributions f± {rj) become similar to 
ribbon-like fettuccine, unchanged along the flow. Since 
/_ results from reverse dynamics, the distributions /-f 
and /_ are elongated along the unstable and stable di- 
rections of the positive time evolution, by the factors e'^*^ 
and e^*^, respectively. Each of distributions /_ and 
consists of two positively and negatively "charged" fet- 
tuccine originating from two oppositely signed separate 
regions of dipole moment distribution /. Since the vector 
field ^ introduced earlier is zero outside the support of 
/, the integration in Eq. (jlSp is restricted to the over- 
lap of three distributions: / and f±. The overlap of 
each of distributions /+ and /_ with / is represented 
by Ni{ti) and Ni{t2) disconnected fettuccine pieces as 
found in Eq. ([T4)) . Therefore the overlap of all three dis- 
tributions is represented by N2{ti,t2) ^ ^^1(^1)^1(^2) ^ 
(Z/a)^e~^^*i"'"*^^ disconnected regions. The linear dimen- 
sions along the isoenergetic shell of each disconnected 
region are estimated as ae~^*^, ae""^*^, and a, which cor- 
responds to the stable, unstable, and fiow directions, re- 
spectively. Therefore the volume of each disconnected 
region is V2{t) ^ a'^e"'^^*^"'"*^^. 

If instead of the second-order response function, we 
were dealing with a three-point correlation function the 
rest of the story would be straightforward. A three-point 
correlation function can be represented in a form of the 
overlap integral /'■^•' (ti, t2) with the differential opera- 
tor £^^dj = {/, •} replaced by the operator of multiply- 
ing by f(r]). In full analogy with the linear response 
case, the typical value of the three-point correlation re- 
sults from cancellations of oppositely signed contribu- 
tions of similar magnitudes and therefore attains a factor 
y/N2V2 ~ g--^(ti+t2)/2 describes the long-time decay- 
ing behavior of the correlation function. The situation 
with the nonlinear response function is apparently more 
complicated since generally the vector field ^ has a com- 
ponent along the stable direction. The derivative of 



a sharp feature in /+ along the stable direction can cre- 
ate exponentially large e^^^ factors. This is the Liouville 
space signature of the exponentially growing components 
of the stability matrix, which affects the response starting 
with second order due to FDT [l^,[l^. The exponentially 
growing components of the stability matrix may seem to 
be a reason for an exponential divergence of the nonlinear 
response, since interaction with the driving field can be 
considered as a kick leading to an infinitesimal deviation 
that grows exponentially with time. This would actually 
happen if the initial distribution was 5-functional con- 
centrated at some point in phase space, and the signal 
was measured by a deviation of the perturbed trajectory 
from its unperturbed counterpart. However, the dipole 
/ that describes the system coupling to the driving field 
is represented by a smooth function. As shown below 
due to this smoothness the exponentially diverging terms 
cancel out completely. This demonstrates that response 
in chaotic systems should be treated using the Liouville 
(distribution-based) representation of classical dynamics, 
while its dual trajectory-based counterpart may lead to 
misleading naive picture. 

To prove the harmlessness of the derivative in the 
second-order response function, we decompose the vec- 
tor field ^ = + ^0 + + into the direction along 
the energy, flow, unstable, and stable components. This 
leads to a natural decomposition of the overlap integral, 
according to Eq. p5)l : 

The term involves a derivative along the stable di- 
rection that provides an additional exponentially small 
e~^^'^ factor, since the distribution is elongated along 

the unstable direction. The dangerous term that involves 

(2) 

derivatives along the sharp feature is represented by /i . 
The aforementioned cancellation can be seen after the 
integration by parts: 

/i") =- j dr7div|_/_/+ - j drjf+eLdjf- . (17) 

The flrst term in Eq. (|17p includes time-independent dis- 
tribution div^_ . The second term contains the derivative 
in the stable direction of the distribution /_ . Therefore, 
it acquires an additional exponentially small factor e~^*^ 
and becomes negligible at long times. Finally, we have 
at long times 

« J df^f_ (-div|_ + eA + ^Edj) /+ . (18) 

All terms in Eq. ([18]) do not have derivatives in the 
stable and unstable directions that provide exponen- 
tially growing or decaying factors. Therefore, according 
to the arguments presented above, the second-order re- 
sponse function has a long-time asymptotic S^'^\ti,t2) ~ 
g-A(ti-i-t2)/2 gijjiiiar to the three-point correlation func- 
tion. Stated differently, the exponentially growing com- 
ponents of the stability matrices that in principle enter 
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the expressions for nonlinear response functions [15|, ll9| , 
do not affect the long-time behavior of the nonlinear re- 
sponse, due to the smoothness of the dipole function /(ry) 
that describes the system coupling to the driving field. 

The qualitative arguments developed in this subsec- 
tions for the lowest-dimensional case n — 1 can be ex- 
tended to a general case n > 2 in a straightforward 
way. This leads to the same asymptotic expressions for 
magnitudes of response functions S^^^(t) ~ e"'*'*/^, and 
S^'^^iti,t2) ~ e"^(*i+*^)/^, where A should be interpreted 
in the sense of Eq. as the sum of positive Lyapunov 
exponents. 



C. Effective negative curvature of configuration 
space 

In subsection IIIIBI we have established a qualitative 
picture of response in classical systems with hyperbolic 
and mixing dynamics. In the forthcoming sections we 
support our qualitative arguments by performing explicit 
analytical calculations of the first- and second-order re- 
sponse functions in a model system of a free particle mov- 
ing along a Riemann surface of constant negative curva- 
ture. Due to strong dynamical symmetry in this system 
(see section |IV] for the details) the response functions can 
be calculated explicitly. In this subsection we present a 
rationale why this simple model system can be viewed as 
a representative example that bears the basic qualitative 
features of chaotic dynamics in a molecule or collection 
of molecules for high enough energies. 

First of all classical motion of a molecular system can 
be represented as classical dynamics of a multidimen- 
sional particle moving in a potential U{r), where r stands 
for a complete set of nuclear coordinates. Trajectories of 
a particle with energy E in arbitrary potential U{r) are 
known to be the same as for a free motion in a curved 
space with the metric = {l-U{r)/E)Sik Al- 
though the potential generates nonuniform motion along 
the trajectories, the dominant feature of the dynamics is 
the exponentially growing separation between the trajec- 
tories in the regions where the metric curvature is neg- 
ative. The regions of the configuration space that cor- 
respond to negative curvature work as defocusing lenses, 
causing instability, i.e. divergence of close trajectories. 
For example, close trajectories diverge every time they 
approach the boundary of the classically inaccessible is- 
land or pass a region of a potential local maximum that 
belongs to the accessible region. Except for occasional 
symmetric integrable cases, passing through the stable 
regions with positive curvature cannot compensate for 
the instability, since stability reflects oscillatory, rather 
than converging features in the dynamics of the trajec- 
tory deviation. In particular, existence of unstable re- 
gions combined with ergodicity ensures exponential di- 
vergence of trajectories over long enough times. 

In addition, when the motion is finite, the accessible 
part of the configuration space at a given energy can be 
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FIG. 2; From the motion in a potential to the motion on com- 
pact surface, (a) For fixed particle energy E, a potential with 
several maxima defines the multiply connected classically al- 
lowed region (b) Idealization of the hard-wall potential: Tra- 
jectories are confined within the original classically accessible 
2D configuration space refiecting on its boundaries. Reflec- 
tions can be viewed as transitions to the antipode surface 
component glued to the original one along the boundaries, 
(c) Trajectories on the deformed smooth version of the re- 
sulting compact surface are shown as solid and dashed lines, 
which correspond to the original and antipode components, 
respectively. 

multiply connected. In the simplest n = 2 case of two 
coordinates the motion occurs inside a disk-like region 
punctured by g forbidden islands (see Fig. [21). Bound- 
aries of the accessible area are represented by curved lines 
where the potential energy U{r) coincides with the total 
energy E. Some fraction of trajectories approaches the 
boundaries so close that this can be qualified as reflec- 
tion. Utilizing the original argument of Sinai [21j reflec- 
tion can be interpreted as continuing motion on the an- 
tipode replica of the accessible region glued to its original 
via the boundaries (see Ref. [H). The resulting compact 
surface has topology of a sphere with g handles (Riemann 
surface of genus g). According to the Gauss-Bonnet the- 
orem, in the 5 > 1 case the average Gaussian curvature is 
negative, which implies regions of instability and results 
in unstable (hyperbolic) dynamics. 

We also note that loose distinction between reflected 
and deflected trajectories that approach the boundaries 
should not be a matter of concern since classical dynam- 
ics constitutes an approximation for a quantum mechan- 
ical problem. The approximation is not valid near the 
boundaries of the classically accessible region, where the 
quantum uncertainty takes over. More formally, one can 
consider a trajectory reflected (or, equivalently, contin- 
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ued on the antipode replica of the accessible region) if 
it penetrates the near-boundary region where quantum 
effects become important. The classical picture of reflec- 
tion becomes exact for the hard wall potential concen- 
trated on the boundaries, e.g. for Sinai billiards [2l[ and 
Lorentz gas. 

Although in a generic system the effective curvature 
depends on position r in configuration space, we will ra- 
tionalize the described above qualitative picture of re- 
sponse in a chaotic system by calculating the linear and 
second-order classical response functions for free motion 
in a Riemann surface Af^ of constant negative (Gaussian) 
curvature. This model that allows an exact solution has 
been serving as a prototype of classical chaos, as well as 
an example of semiclassical quantization [ij, [ij, [22| • 

IV. GEOMETRY AND DYNAMICS ON A 
RIEMANN SURFACE WITH CONSTANT 
CURVATURE 

A. Free particle on a Riemann surface 

In this section we describe classical dynamics of a free 
particle moving along a Riemann surface of constant 
negative curvature (Gaussian curvature). This will be 
done by making use of strong dynamical symmetry. This 
system is one of the best studied models of chaotic dy- 
namics. Although detailed reviews of its symmetry exist 
in the literature (see e.g. Ref. |23 ). we summarize the 
necessary results for the sake of completeness and fur- 
ther extension to the case of Langevin dynamics. Noisy 
dynamics associated with geodesic flows on surfaces with 
constant negative curvature is considered in detail else- 
where Hi]. 

The Lagrangian of a free particle contains only kinetic 
energy: L = mgikr^r^ /2. The corresponding unper- 
turbed classical Hamiltonian 

= ^^^^ 

does not depend on x if expressed in terms of the ab- 
solute value of the momentum C,. Here we introduced 
covariant and contravariant metric tensors gik and g^^ . 
The curvature (also referred to as the Gaussian curva- 
ture) of configuration space is expressed in terms of 
derivatives of the metric tensor. 

Since unperturbed classical dynamics conserves energy, 
the value of C is preserved, and evolution actually oc- 
curs in the reduced phase space (energy shell) . If 
we are not interested in a trivial case of zero energy, dis- 
tributions can be thought of as functions p{x,^), with 
X G and C > is the absolute value of the particle 
momentum. A compact 3-dimensional smooth manifold 
is the subspace of the phase space that consists of 
points with unit length of the momentum vector. Points 
X G are specified by two coordinates r and momen- 
tum direction angle 0. The dynamical symmetry for the 



free-particle dynamics originates from the fact that there 
is a smooth action of the group G = 50(2, 1) in the 
reduced phase space of the system, consistent with the 
dynamics. The details are presented in Appendix[Xl The 
dynamical symmetry has a simple and clear interpreta- 
tion in infinitesimal terms, i.e. action of the correspond- 
ing Lie algebra so(2, 1), where the algebra generators are 
implemented as vector fields in M^. 

As discussed in Section |lTl the evolution in the phase 
space is determined by equation r) — Lrj. In this equa- 
tion Lt] is a vector field that belongs to the tangent space 
of phase space point rj. Hereafter we adopt an agreement 
used in differential geometry by identifying a first-order 
differential operator of differentiating along the vector 
field with the vector field itself. Due to energy conserva- 
tion, vector field L is tangent to the reduced phase space 
M^. We introduce a vector field cti that determines the 
geodesic flow and generates the translation along trajec- 
tories, so that the Liouville operator of free motion is 

L = Cai. 

Another natural vector field in is CTz which corre- 
sponds to the generator of the momentum rotation and 
keeps the coordinates unchanged. This operator can be 
represented as <7z — d/dO in terms of a derivative with 
respect to the particle momentum direction 9. 

Operators ci and (Jz can be easily expressed in terms 
of r and p. Performing the variables transformation we 
obtain: 

dpi dvi dn dpi dpi 

Here Eij is the asymmetric tensor, £^12 = — -E21 = 
(det g**"')^^/^, £'11 = E21 — 0. Both operators cti and 
CTz obviously conserve energy: uiC = o-zC, — 0. 

Finally, we introduce the third vector field in Al^ as 
o'2 = [ci , cTz] where the commutator of vector fields is 
understood as the commutator of the corresponding dif- 
ferential operators. A simple calculation shows that in 
the case of constant negative curvature the vector fields 
CTz, (Ti, and (72 form the Lie algebra so(2, 1) with respect 
to the vector field commutator (see Appendix [B] for more 
details): 

[ci, 0-2] = (7z , [0-1,0-2] = 02 , [ct2, CTz] = -O-l . (20) 

The group action of SO {2, 1) in the reduced phase space 
is obtained by integration of the algebra action. 
Any compact surface of constant negative curvature 
may be represented as a result of factorization of the 
hyperbolic plane with respect to translations that con- 
stitute a finitely generated infinite discrete subgroup of 
50(2, 1) (see Appendix E]) . The entire hyperbolic plane 
possesses the complete SO {2, 1) symmetry, which makes 
the infinite motion there integrable due to the presence of 
the integral of motion similar to the angular momentum 
in the case of 50(3). The factorization destroys global 
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SO{2, 1) symmetry and folds back the trajectories, which 
renders the dynamics strongly chaotic. 

The Poisson bracket can be computed in a standard 
way, making use of the fact that locally it coincides with 
the Poisson bracket for a free particle moving in the en- 
tire hyperbolic plane H (see Appendix [BJ . The canon- 
ical Poisson bracket may be rewritten in terms of the 
linear differential operators introduced above acting on 
two functions f{x,() and g{x,() as 



{/,5} = 



(21) 



The Poisson bracket in Eq. (|2ip is expressed as a bilinear 
form of generators of a Lie algebra with constant coeffi- 
cients. This is is another important manifestation of the 
dynamical symmetry in the problem of free motion on 
the surface of constant curvature. 

The form of the Poisson bracket as well as the com- 
mutation relations suggest that the Lyapunov exponent 
is equal to A = \J —K. This can be clearly seen from 
the Jacobi equation y -\- Ky = for the magnitude of 
the normal component of the deviation from a given 
geodesic [11]. The stable and unstable directions are 
given by linear combinations of and a^. The dis- 
placements along the flow ai are conserved. The de- 
viation of the trajectory at < = can be written as 
S'qiQ) — a{q), where cr is a first order differential op- 
erator corresponding to a certain direction in the phase 
space. The evolution of the phase point 77-1-1^77(0) is given 



by e^^^T] + arf) = 'q{t) 



ri{t). We see that 



after time t the deviation from the phase point 77(i) is 
determined by the vector field a{t) = e°'i*cre~'^i'. Com- 
mutation relations (I20|) allow to find this Heisenberg rep- 
resentation of any operator decomposed over the basis 
elements of so(2, 1). In particular, we find that 



-o-it _ „±t 



((72 ± (7z) : 



(22) 



and thus conclude that the local stable and unstable di- 
rections are given by (T2 — (Jz and 172 + tXz, respectively. 
The fact that the form of these vector fields is conserved 
by the dynamics is an important manifestation of the 
dynamical symmetry. 

Dynamical symmetry also leads to the following gen- 
eral relation: 



/ dxaif{x) = 

JM3 



(23) 



for Z = 1, 2 or 2 that is valid for any smooth function 
f{x) where dx in is an invariant integration measure 
with respect to the SO{2, 1) action (see Appendix 1X1 for 
details). Eq. (j23[) implies an integration- by-part rule 

dxf{x)aig{x) = - dx {cjif [x)) g{x) ., (24) 

that will be an important ingredient of our analytical 
calculations. 



In conclusion we emphasize that the dynamical sym- 
metry with respect to the action of the group G = 
S0{2, 1) does not mean symmetry in a usual sense, i.e. 
that the system dynamics commutes with the group ac- 
tion, but rather reflects the fact that the vector field L 
that determines the classical dynamics is represented by 
an element of the corresponding Lie algebra so(2, 1). This 
allows to apply representation theory as described below. 



B. Decomposition of tlie free- particle dynamics 
using irreducible representations 

The smooth action of G in can be interpreted 
as that the space H of smooth functions in con- 
stitutes a representation of G, which turns out to be 
a unitary representation (see Refs. [H, [H, [H and Ap- 
pendix[X|, and therefore can be decomposed into a direct 
sum of irreducible representations of G. The spectrum 
SpeCg(Af^) C G of the Riemann surface is defined as a 
discrete subset of the space G of irreducible unitary rep- 
resentations of the principal series that participate in the 
decomposition of functions in into irreducible repre- 
sentations: 



(25) 



seSpcco(M2) 



The spectrum SpecQ(M^) consists of imaginary numbers 
s that characterize representations of the principal series. 
The effects of the inclusion of other (complementary and 
discrete) series are discussed after their review in Appen- 
dices El and The evolution in the space of distribu- 
tions is decomposed into a set of uncoupled evolutions 
that correspond to relevant irreducible representations. 
The component dynamics is determined by the reduced 
Liouville operators: 



dgs{x,C;t) 
dt 



= ~'CL{s)gs{x,C,t) . 



(26) 



where L(s) corresponds to an irreducible representation 
labeled by s. The distribution components can be further 
decomposed as 



(27) 



k— — c 



using the eigenstates of the angular momentum operator 
az- These satisfy the following properties: 



azipk{x;s) = iktljk{x;s) 



(28) 



a±iJkix; s) = [±k + -~ s] 'ipk±iix; s) , 



where we introduced the raising and lowering operators 
(7± ~ ai ± 1(72 that are anti-Hermitian conjugated, i.e. 
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The description of unitary representations of G, re- 
viewed in detail in Appendix [Ul (see also Refs. [2^. l25h. 
allows to represent the eigenstates ipkix; s) by the func- 
tions on the circle ^fc(u) for any given s € Spec(M^). In 
the case of imaginary s (principal series), when \E',„(m) 
constitutes an orthogonal normalized set with the nat- 
ural scalar product, this leads to normalized functions 
ipk{x;s). Identification of the functions ijjkix) with 
^'fc(u) establishes an isomorphism between two irre- 
ducible representations, the first being Ti^ that partic- 
ipates in the decomposition of Eq. (jA4p , the second be- 
ing its standard representation in function in a circle de- 
scribed in Ap pendix [Cl According to the Shur lemma 
(see, e.g. [23]) the identification (isomorphism) is deter- 
mined up to a factor. Its absolute value can be fixed 
by requiring that the function 'ipo{x;s)^ identified with 
^^oiu) = 1 is normalized. To fix also its phase we note 
that ipoi^', s) does not depend on 9 and turns out to be 
an eigenfunction of the Laplacian operator in 
as described by Eq. (|C5|) . Since the Laplacian is a real 
operator, the eigenfunction ipoi^] s) can be chosen to be 
real. Hereafter, we implement an agreement that '>po{x; s) 
identified with ^o{u) = 1 is real. This determines the 
functions up to a sign, the latter being of no importance: 
*fc(u) = e''^". 

The form of ai in the angular representation is fixed 
up to a phase [24] : 



L{s) = (Ti = sin u 



1 - 2s 



du 



(29) 



Thus, the original dynamics of a distribution gs{x;s) in 
Eq. (P7)) is mapped onto an effective classical dynamics of 
a distribution Giu) = X]fc^-oo Ps,fc(C)^A;(u; s) defined on 
a circle with the Liouville operator (|29p . The resulting 
effective problem is one-dimensional and can be solved 
exactly. 

The term that describe the interaction with the driving 
field can be also decomposed in irreducible representa- 
tions. The coupling (polarization) / is represented by a 
function of the particle position only (a function in Af^) 
or, stated differently, a phase-space function that does 
not depend on the particle momentum. It can be equiv- 
alently interpreted as a function in independent of 6, 
hence a^f = 0. Thus, / can be expanded as 



/ 



seSpcC(,(M2) 



Bsi^oix; s) 



(30) 



The sum in Eq. (j30p runs over the spectrum of the com- 
pact surface corresponding to the principal series repre- 
sentations characterized by imaginary s (see Appendix 
[C|) . We emphasize that, as noted earlier, the functions 
ipoi^', s), being actually functions of the particle position 
only, are the eigenfunctions of the Laplacian operator 
on the Riemann surface. Typically the dipole / is 
a slow function of coordinates, and therefore only few 
terms provide substantial contributions to the expansion 
of Eq. ([50)1 . We also note that, although the Laplacian di- 



agonalization on an arbitrary Riemann surface with con- 
stant curvature is a complex problem, its eigenfunctions 
are tangible and intuitively simple objects. Combined 
with the previous results, the expansion ([50]) leads to 
closed analytical expressions for the response functions 
of the original problem. 



C. Effective dynamics on the circle 

We have identified the dynamical symmetry that al- 
lows to map the original chaotic dynamics on onto a 
tractable dynamical problem on the circle. In this subsec- 
tion we consider reduced dynamics that corresponds to a 
principal series representation labeled by s with Ims > 0. 
In the following calculations of the response functions, we 
will need the expansion of e~^*^o{x; s) over basis vectors 



-Lt 



ipoix; 



^ Ak{t;s)'4;k{x;s) 



(31) 



fc=-c 



In particular, the linear response is expressed in terms 
of the coefficient Aq. These coefficients that are actually 
matrix elements of the evolution operator between the 
angular harmonics have been calculated earlier [l^ . 
The result immediately follows from the description of 
irreducible representations using a construction of an in- 
duced representation (see, e.g. Ref. IH). The derivation 
presented below allows an extension to the Langevin dy- 
namics, associated with the original classical dynamical 
problem. The picture based on the Fokker-Planck equa- 
tion will be presented in detail elsewhere [l^ . 

The coefficients Ak{t; s) are calculated by implement- 
ing the representation on the circle, introduced above. 
Since e~^'^'*-'*^'o(u) represents a function obtained as 
the result of the evolution operator action on '^o{u) = 
1, it can be found by solving the dynamical equation 
dtg{t,u) + Lg{t,u) — 0. For the principal series repre- 
sentation labeled by s, the equation adopts the following 
form: 



dtg{t,u) 



svau 



d_ 

du 



1 - 2s 



cosu g{t,u) = 0. (32) 



The solution of the first order partial differential equation 
supplemented with an initial condition g{Q, u) = 1 can be 
found using the method of characteristics: 



-L(s)t.- 



e "^■■''^5'o(u) = e 



= (cosh t + sinh t cos u) " 



u 




cos- 









1 



-2t, 



(33) 



At long times t > Q the solution is concentrated near the 
stable stationary point u = tt. This reflects the collapse 
of the reduced phase space distribution function along 
the stable direction. Indeed, the width of the region 
where e~''^^'*^*^'o(u) is not exponentially small vanishes 
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as cx e~*, according to the fact that for the chaotic dy- 
namics on Af ^ with Gaussian curvature K the Lyapunov 
exponent is equal to \J ~K. 

The expansion coefficients in Eq. (|3ip are given by 

^fc(a;s)= y'^**(u)e-C^(^)*vl/„(u). (34) 

The integral can be reduced to a standard integral 
representation of the Gauss hypergeometric function 
2F\(a, 6, c, z) defined as a series 

,1^ r(a)r(6)r(c + n)n! 

r(z) being the gamma function [26,] . An exact expression 
for fc > reads: 



U + 77 - 



.r(fc 



fc!r(i-.) 

s, fc 



, t 

tanh — 
2 



(36) 



1,-sinh^ - 
2 



Various representations of these coefficients using other 
special functions have been derived in the context of two- 
point correlations [isl. [2^. 

Since solution g is even, g(i, — u) = g(t^ u), coefficients 
Ak are symmetric, A^k = Ak- The time reversal prop- 
erty 



g[-t, u) ^g(t,u + Tr). 



(37) 



that follows directly from Eq. ()32p . immediately implies 
Ak{—t) = {—!)'' Ak{t). We can represent Ak{t;s) in an 
alternative form suitable for studying its long-time be- 
havior 1261: 



Ak{t;s) 



2(-l)^-(l-e-^^)''r(fc+l-.) ^_,,, 
V^r(i - s) 



X (38) 



Rc 



r(s)e^ 



r(fc 



2i^i(fc + i-s,fc+i,l-s,e-2*jy 



Its expansion in e^-^* using Eq. ((35)) corresponds to the 
spectral decomposition of the regularized evolution oper- 
ator for the irreducible representation labeled by s. The 
properly interpreted eigenvalues of the Liouville opera- 
tor are known as Ruelle-PoUicott resonances for chaotic 
systems. They can be obtained using a physical ap- 
proach. Adding weak Langevin noise results in adding 
an infinitesimal second-order differential diffusion oper- 
ator to the classical Liouville operator. This results in 
a Fokker-Planck operator £ = — kV^ -I- L that repre- 
sents a regularized version of the Liouville operator [27| . 
The resonances are more often treated as mathematical 
objects in a form of generalized eigenfunctions in a prop- 
erly chosen rigged Hilbert space (see e.g. Ref. [isl ). The 
physical approach to spectral decomposition of linear and 
nonlinear response functions for a free particle moving on 



a surface of constant negative curvature has been devel- 
oped in Ref. [28l . 

We will be interested in long-time behavior of the re- 
sponse functions and expect that the asymptotics 



Akit;s) 



2(-l)'=r(fc + 



V^ra - s) 



-^e-*/2Re ^ 



\r{k 



s) 
(39) 



at i — !■ CX) are relevant. However, we will see that, 
since the approximation breaks down for higher harmon- 
ics fc ^ e* and does not vanish at fc — > cxd, it cannot be 
immediately utilized for nonlinear response calculations 
where the result is given by the infinite series over fc. In 
the region t > 1 and 1 <C fc ^ e^* one can use an approx- 
imation for Ak{t) in terms of the MacDonald function 
(modified Bessel function of the second kind): 



Ak{t;s) 



2(-l)'^ 

v^r(i-,s) 



e-'^^k-''Ks{2ke-'). (40) 



V. CALCULATION OF RESPONSE 
FUNCTIONS 

In this section we substantiate the qualitative picture 
of response presented in Section IIIII which is based on 
general expressions ([TOl [TT|) for a particular model of 
free motion on a compact surface with constant nega- 
tive curvature. Heretofore we have not made any as- 
sumptions about the shape of the initial distribution po 
in the general expressions for response functions. In the 
case of free motion, the equilibrium phase space density 
Po depends on the momentum absolute value only. The 
fluctuation-dissipation theorem (FDT) relates the linear 
response function to the correlation function in a system 
at thermal equilibrium. To reduce bulky calculations, 
in both linear and second-order response functions, we 
employ an analog of FDT and notice that the action of 
the evolution operator on the most inner Poisson bracket 
can be reduced to the time derivative. It allows recasting 
Eqs. (|10lll|) in a form: 



5(i)(ii)=9i, / drjf{ri)e-^'^f{ri) 



dpo 

dc 



(41) 



(42) 

The general expressions given by Eqs. (j4ip and (|42l) ap- 
ply to any equilibrium distribution, including the canon- 
ical Po cx e"^^, microcanonical po cx S{E — Eq), and 
cumulative microcanonical OePo oc S{E — Eq) distribu- 
tions. To obtain the final expressions in the most trans- 
parent form we will use the cumulative microcanonical 
distribution with the cut-off energy set to Eq = 1/2. 

We are able to approach an apparently intractable 
problem of calculating the response in a chaotic sys- 
tem because of strong dynamical symmetry present. 
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The preliminary steps of our calculation have been de- 
scribed in the preceding section. Dynamical symmetry 
leads to decomposition of the original dynamical prob- 
lem into a set of uncoupled evolutions that correspond 
to irreducible representations of the dynamical symme- 
try group 50(2,1). The dipole moment / is expanded 
over momentum-independent basis functions 4'o{x; s) ac- 
cording to Eq. (|30p. where s labels relevant irreducible 
representations. The relevant matrix elements Ak(t;s) 
of the evolution (Perron-Frobenius) operator are given 
by Eq. The coefficients Ak{t;s) constitute the dy- 

namical part in the calculation of the response functions. 
The geometrical part of the calculation is represented 
by the integrals over the reduced phase space, hereafter 
referred to as (geometrical) matrix elements. The ge- 
ometrical part turns out to be trivial in the linear re- 
sponse case while the calculations for the second-order 
response are more involved. Computation of the matrix 
elements involved in the nonlinear response is one of the 
main technical results of this manuscript. The second- 
order response function will be obtained in a form of an 
infinite series. A nontrivial character of the series con- 
vergence reflects a nontrivial way how the exponentially 
growing terms cancel out in the nonlinear response, as 
described in Section UlI Bl 

In this section for the sake of simplicity the calcula- 
tions are performed in the case when the dipole moment 
/ has a single component / = ipaix; s) in the expansion 
over irreducible representations ([50)1 . Stated differently, 
the dipole function / is represented by a single eigen- 
mode of the Laplace operator. Some details for the case 
of two components, / = Bs^il^oix; si) + Bs^ipoix; S2) (su- 
perposition of two Laplacian eigenmodes), are presented 
in Appendix [Dl Numerical results for the second-order 
response function in this simplest case that involves dif- 
ferent resonances are given in Section IVII 

A. Linear response 

Since the basis functions ^k{x]s) with different k are 
mutually orthogonal, the linear response function is de- 
termined by the dynamical part alone and can be repre- 
sented in a form 

00 

5«(t) = |/dCAo(CM)§^, (43) 


with 

A(i;s) = -^e-*/2x (44) 

according to Eq. ((38l) . For large t the linear response 
function shows oscillatory decay as e"^^/^*"'*. The ex- 
pansion in powers of e~^* can be interpreted as the expan- 
sion over Ruelle-PoUicott resonances. The resonances can 



be found as the eigenvalues of the regularized Liouville 
operator [2^]. Only even resonances u}2k = 2k + 1/2 ± s 
contribute to the response function. The expansion is a 
convergent series over e"^*^ if e~^*^ < 1, i.e. ti > 0. 

In a general case the coupling f{x) can be decomposed 
in irreducible representations labeled by s [ see Eq. (ISH]) ] , 
and the linear response function becomes a linear com- 
bination of contributions with the coefficients B^. 
Due to orthogonality of the basis functions, there is no 
coupling between different representations. 

B. Second-order response: matrix elements 

The second-order response function for a free parti- 
cle on a Riemann surface with constant negative cur- 
vature is obtained by substituting the specific expression 
for the Poisson bracket [Eq. (HH)] into the generic expres- 
sion given by Eq. ([1^ . We further simplify the calcula- 
tion by making use of the evolution operator unitarity 

^g-it2^ — git2 -vvith respect to the natural scalar de- 
fined by Eq. (|C2p . Performing the integration over by 
parts and making use of 9po/i9Clc=o,oo = results in: 

S'-^Hh,t2) = -g^J dxdCx (45) 

(^2 (ale'^^^*^/)*(al/)(e--^';*V) 

-|-(e^^«*V)*(^i/)(e-"^^*^/) 

+ (e-'^*V)*(a2/)(^.e--';*V))^. 

The transformation to Eq. (|45p is an important step 
towards expressing the result in terms of the quantities 
calculated in the previous sections. This allows to avoid 
propagating the vector fields (T2 and az, and instead deal 
only with the action of the evolution operators e^*^ and 
e~^*i on the dipole moment /. The latter is decomposed 
according to Eq. (j30p in zero harmonics ipoix] s) of irre- 
ducible representations labeled by s, i.e. the Laplacian 
eigenmodes. Then we use the expansion pip to express 
e~^*ipo{x; s) in terms of the angular harmonics ipk{x; s) 
which constitutes a basis set in the representation. 

Now we approach the most challenging part of the cal- 
culation represented by the integration over the reduced 
phase space. Generally, according to the decomposition 
of /, the integration may involve functions from different 
irreducible representations. In this section we consider 
the simplest case / = ipoix; s) of the dipole function rep- 
resented by a single eigenmode of the Laplacian. A more 
complicated case of two representations is treated in Ap- 
pendix |D] 

One needs to find a way to calculate the integrals of 
triple products like ip^{x; s){ai^24'o{Xy s))tpmix; s) over 
the reduced phase space M^. This matrix element in the 
second-order response function is essentially less trivial 
part than its counterpart in the linear response function. 
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where the integral of the double product of 4'n{x; s) has 
been done based on orthogonality and normalization of 
the functions. The problem arises from the fact that in- 
tegrals of relevant triple products over the reduced phase 
space do not have natural representation in terms of the 
effective system on the circle. Nevertheless, we will show 
that this geometrical part of the second-order response 
can be worked out using relatively simple tools, and can 
be expressed in terms of few quantities related to the 
Laplacian eigenmodes. 

The integration in reduced phase space x — {r,6) is 
performed over the particle coordinates r e and the 
locally defined angle 9 that represents the 2D momentum 
direction. The local (for fixed r) dependence of ijjn{x; s) 
on 9 is given by ipnix^s) oc e™^ due to the definition 
dgipnix; s) = in^n{x; s) of the basis functions. The op- 
erators iJi and a2 applied to ijjQ{x;s) in the middle of 
the integrand are expressed via the ladder operators a± . 
Therefore, integration over 9 in Eq. (|45p results in zero 
unless the integrand is locally independent of 9, i.e. 



between the matrix elements in an explicit form: 



Af3 



dx {ipnix; s))* {cr±il;o{x; s))iprn{x; s) cx (5„,m±i .(46) 



We now introduce the matrix element Ofe and bk as 



Ofc = / dxiplix; s)'ipa{x; s)ipk{x; s) , (47) 
bk=2 dxipl{x;s){a+^oix;s))^k-iix;s). (48) 



As established in Section IIVI the differential operators 
a± can be moved from one part of the integrand to an- 
other one, and thus one can integrate by parts [Eq. ([24]) ]. 
Integrating by parts in Eq. (j48p and making use of the 
properties (|28p of cr operators we derive the following im- 
plicit recurrence relations: 

(49) 

2 J dx {a^ipk+iyipntpk - 2 / dx ^pl^^^^pocr+^pk 
{2k + I - 2s){ak - ak+i) , 
and 



(50) 



bk+i =2 dxil)l^-^{(j+il)Q)ipk 



bk+i ^2 1 dx ■iljl_^_j^{a+'ipo)ipk 
4 



2k + l + 2s 

4 



2k + l + 2s 



dx (o-+-0A;)*(o'+'0o)V'fe = 

dx'iljl{{a^a+i}jo)'ipk + ia+tpo)<J-'4'k) 



l-4s2 2fc-l + 2s, 

-«fc + ^, , , , ^ Ofc ■ 



2fc -t- l-f 2s 2k + l + 2s 

We further make use of Eq. ([49]) to express ak+i via 
and bk+i, and applies Eq. (jSOp to find bk+i in terms of 
Qk and bk- This allows to recast the recurrence relations 



4:k{k+l) 



2fc - 1 + 2s , , , 
rflfe - TTTT— rT^5 —bk, (51) 



(2fc-hl)2-4s2 " (2/s 1)2 - 4s2 
l-4s2 2fc-l-^2s, 

-Ofc + Ofc . 



2k + l + 2s 



2k + l + 2s 



The matrix elements bk can be excluded to connect three 
adjacent matrix elements afe_i, and ak+i'- 



8/fc2 



l-4s2 



(2fc + l)2-4s: 



{2k_-rf_^ 
(2fc + 1)2 -4s= 



ak-i . (52) 



As mentioned earlier ipQ can be chosen real, which leads 
to the relations ai = a_i = ao/2 (see Appendix [Dl for 
the details). This allows to solve the recurrence relations 
and express all matrix elements in terms of a single real 
number aq. In addition, in Appendix ID] we consider a 
more complicated case of matrix elements composed of 
functions that involve different representations. 

Putting all contributions to the second-order response 
function together, we can write it down in the following 
form: 



||f;(-l)"(a„-a„+i)(n+i + s) x 

^ (<2 (C^i ) a: (^2) - K iCh ) A„+i (Ct2 )) ) 

\n + l)A^+,iCh)AUCt2) + nA;(Cti)A„+i(Ct2) 

(53) 

Here we used the short notation An{t) for A„(i;s) and 
the fact that the product {2n+l — 2s)A'!^^i{t2; s)An{ti; s) 
is real. We used the symmetries A_„(i) — An{t) and 
a_„ = a„ discussed in Section flV CI and Appendix [Dl to 
restrict the consideration to non-negative n. 

We note that the coefficients A*(— t2;s) at negative 
time — ^2 < enter the expansion for the backward evo- 
lution of ipoix', s), which appears to be crucial for the cal- 
culation below. The time reversal symmetry A„{—t; s) = 
{—l)^"An(t;s) was also used in Eq. 

The result for the second-order response function ([53]) 
is written in the form of an ordinary series of terms whose 
expressions are either completely known analytically (see 
Eq. (|36p ) or recurrently expressed via a single number 
ao that characterizes function ip(){x;s). The functions 
^po do not depend on the angle 9, i.e. ipo{x) = (j){r;X), 
where, according to Eq. (jC6p (/)(r; A) is an eigenmode of 
the Laplace operator V2 with the eigenvalue A on our 
Riemann surface. Therefore, Eqs. (|A3p and (jC6p allow 
to express the only undetermined parameters s and ao 
in Eq. ([55|l in terms of the relevant eigenmode (/)(r; A) of 
the Laplacian: 



ao 



M2 



dr (0(r;A))^ 



iV-4A- 1 



(54) 
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Finding the eigenmodes of the Laplacian analytically is 
a complex problem. Fortunately this is not necessary for 
our purposes. It is known that for a given constant neg- 
ative curvature K < there is a family of non-equivalent 
Riemann surfaces of genus g that can be parametrized by 
a space Aig with the dimension dimA^g = 4g — 2 known 
as the moduli space. In particular, ap and s can be viewed 
as some nontrivial functions on the moduli spaces Mg, 
and computation of these functions is a complex prob- 
lem. Yet, we can note that for any given values of ap 
and s we can find a Riemann surface that implements 
them according to Eq. ([5^ . Therefore, we can treat oq 
and s as free parameters that characterize some Riemann 
surface of constant negative curvature. 

In a more general case, when the dipole / is represented 
by a finite superposition of Nf Laplacian eigenmodes, 
the second-order response can be expressed in terms of 
a larger set of parameters sj and Cq'^'*'', with i,j,k — 
1, . . . ,Nf as described in Appendix ID] We can apply the 
same argument to treat them as free parameters for a 
Riemann surface of a high enough genus g. 

C. Second-order response: summation of the series 

We have expressed the second-order response in terms 
of an ordinary converging series for the cases of a sin- 
gle Laplacian mode Nf = 1 [Eq. ((53)) ] and two modes 
Nf = 2 [Eq. (|D7p ] , respectively, in the expansion of the 
dipole ([50)1 . A general expression for a finite number Nf 
of modes has a similar structure. Each term of the series 
is known analytically in the form that involves special 
functions and solutions of recurrence relations. Or goals 
are to derive the long-time asymptotics of the second- 
order response S^^\ and develop a computationally effi- 
cient procedure for S*^^-* (ii, ^2) at finite times. 

We start with the asymptotic behavior of the second- 
order response function for large ti and t2. Naively, one 
would plug the long-time asymptotic of A„(i; s) from Eq. 
(15^1) into Eq. However, this asymptotic of An{t;s) 

is valid only for fixed n, if ne~* <C 1. This can be con- 
firmed by comparing the consecutive terms in the expan- 
sion of the hypergeometric function in Eq. (j36p in powers 
of e~*. Moreover, although An{t; s) represent the Fourier 
expansion coefficients of a smooth function p3p . their 
asymptotic ([M)) does seem to vanish as n ^ 00. Count- 
ing the powers of n in the second part of the summand 
in Eq. (|53p estimated using these asymptotic expressions 
reveals that the resulting series fails to converge. 

Indeed, the series (j53p can be represented in the form 

oc 

^^'^(^1,^2) = Y.{-irF{n;h,h) , (55) 

n=0 

where the dependence of F{n; ti, t2) on n is slow for n ^ 
1. We may find the long-time behavior of the terms as 

oc 

F{n;h,t2) = J2 Fki{n-MM)e-^'''''^''' ■ (56) 

kl=0 



The double expansion originates from the decomposition 
of An{t]s) in powers of e~^*. The first term in the lat- 
ter decomposition, specified by Eq. turns out to be 
independent of n in the limit of large n, apart from irrel- 
evant slow quasi-oscillations oc that cannot affect the 
series convergence. A naive long-time asymptotic of the 
series is given by ^jj(— l)"i^oo(?i; ^i, ^2), where the n-th 
term is estimated as 

Foo(n;ti,i2) = »^ (fln - a«+i)e ^ 
X Re(C„(s)e''*i)Re(i:'„(s)e''*^). (57) 

where quasi-oscillatory dependence of a„, Cn{s) and 
Dn{s) on ri ^ 1 does not affect the convergence. Taking 
into account the asymptotic form of a„ [see Eq. (jEip j. 
we conclude that the series is obviously divergent as 
E„(-l)"^oo(";ti,i2) ~ E„(-l)"v^- In what follows 
we show how to overcome the apparent divergence and 
calculate the asymptotic of the second-order response 
function analytically. 

The series in Eq. ([55]) that represents the response 
function must converge for fixed values of ii and t2. 
Indeed, the terms F{n;ti,t2) oc exp[— 2n(e^*i -f e^*^)] 
vanish exponentially, although only for extremely large 
n ^ e*i,e*^ ^ 1. The decay rate is determined by the 
fact that g{t,u) = e~^*'i>o{u) in Eq. ([55]) is smooth on 
scales that do not exceed the smallest fettuccine size e~*. 
The ultimate convergence allows to safely regroup the 
terms of the series: 

00 

^(-l)"F(n) = -F(0)+ (58) 

1 

- ^(F(2fc) - 2F(2fc -f' 1) -I- F{2k + 2)) . 

fc=0 

After the terms are regrouped the initial naive ap- 
proach based on approximating F{n;ti,t2) by its long- 
time asymptotic -F'oo('^; ^ii ^2) results in a converging se- 
ries, since the linear combination in the summand rep- 
resents the discrete counterpart of the second derivative 
(P F^(i{n) / d-n? which decays as oc n~^l'^ . According to 
Eq. ([57)1 all terms in the resulting series have the same 
time dependence. Therefore, the long-time asymptotic of 
the response function S^'^^(t\,t2) is represented by a su- 
perposition of four components e~(*i+*^)/^='=*(*i='=*^' with 
the coefficients in the form of convergent series. 

Coefficients Fki(n; ti,t2) at higher orders in the expan- 
sion (I56p grow faster with increasing n as Fki{n; ti, 12) ^ 
n^("+*'''Foo(n; ti, t2)- The series re-grouping approach of 
Eq. ([58)1 can be generalized to eliminate the apparent di- 
vergences for higher-order terms in Eq. (|56p . and obtain 
the higher-order terms in the asymptotic expansion of 
the second-order response in powers of e~^*^'^. Instead 
of that we suggest an alternative procedure that allows 
(i) to demonstrate the existence of a long-time asymp- 
totic expansion of S''-^-' in powers of e"^'^ and e~^*^, (ii) 
derive relatively simple expressions for the expansion co- 
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efficients in any order, and (iii) develop an efficient nu- 
merical scheme for computation of the response at finite 
times. 

We start with deriving the asymptotic expansion. To 
that end, provided ii, i2 S> 1, we introduce an intermedi- 
ate ^ 1, so that for n > N the terms F{n;ti,t2) 
are represented by smooth functions of n. We fur- 
ther partition the sum S of the series into the finite 
sum Sn = Sn<jv(~l)"'^('^) ^^'^ remainder i?jv = 
l)"F(n), followed b y ev aluating the remainder. 
Implementing the definition [26| of the Eulcr polynomi- 
als En{y), we derive the following identity, based on the 



Taylor expansion of (1 -I- e 



m X 



d/dz: 



F{z+l)-Fiz)^J2 



2n! dz" 



(F(z + 1)-F(z-1)), 



which is used to calculate the sum of consecutive terms 
pairwise, so that the remainder Rn = X]n>Ar(~l)"^("-) 
of the almost alternating series becomes related to the 
first term F{N) that is not included in Rn, and its 
derivatives: 



R 



m=0 

2 ml \ 



1 



i-l)^+^{-F{N) + -F'iN) + 



(59) 



We further introduce the following "improved" partial 
sums: 



"at — JN 



^-^^f:^F(-HN), (60) 



2 ^ — ' m, 

m=0 



so that 5(2) = pf^ + 0(F(*^+i)(7V)). Due to the 
smoothness of F{n) the deviation of 5*^^^ from its im- 
proved approximation P^*^'* may be estimated as ^ 

F(7V)/iV(*^+i). Since P^^^ are determined by F{z) for 
z < N , & choice of <C e*^, e*^ allows to use the expan- 

(M) 



sion of Eq. (|56|) which leads to an expansion of P, 



N 



in powers of e and e with 



N 



p: 



(M) 



Y.{~itF^,^,{k-tiMy 



(62) 



fe=0 



2 ^ mi 

m=0 



The reason for introducing the improved partial sums 
P^^^ is that the expressions for the expansion coefficients 



(M) 



have finite limits at A/' — > oo provided M > 



p. 

mi+m2- Stated differently, the second term in Eq. ([62]) 
can be viewed as a set of counter-terms that eliminate 
the divergence in the first term. In the limit N ^ oo 
Eq. leads to the expansion 

S^'Htut2)= V 5^^U(^i'*2)e-^^'"^*^+™^*^) (63) 



mim2 



with 



s^^U{h,t2) = \im p^r:x:::\h,t2) 



(64) 



The expansion of Eq. (|63p is represented by an asymp- 
totic rather than converging series. It can be viewed 
as a double expansion of the second-order response func- 
tion 5'(2)(ti, i2) in Ruelle-PoUicott resonances which orig- 
inates from the uncoupled evolution of the chaotic sys- 
tem during time intervals ti and ^2- This will be demon- 
strated explicitly [28] by analyzing the noiseless limit of 
the corresponding Langevin dynamics. 

At this point we should note that our derivation of the 
asymptotic expansion has been somewhat frivolous. First 
of all, since the terms F{n; ii, t2) of our original series in- 
volve the coefficients am , we do not actually know a func- 
tion F{z; ti, t2) that represents the series, only its values 
for positive integer integer z — n are available. In partic- 
ular the derivatives d^Fmim2{N) in Eq. ([62)1 have been 
not defined yet. Second, in obtaining the asymptotic 
expansion we were not controlling the neglected terms. 
This is especially dangerous in our case since we derive 
the complete asymptotic, which involves computing the 
terms that are exponentially small compared to more se- 
nior terms in the expansion. We need to make sure that 
the terms that are neglected in computing a certain ex- 
pansion coefficient are small also compared to the higher 
terms that are kept in the asymptotic expansion. 

A sketch of the appropriate derivation is presented in 
Appendix [E1 In particular, we demonstrate that the co- 
efficients a„ that enter the expressions for F{n;ti,t2) 
can be represented by some analytical function a(z), 
and therefore the terms F{n; ti, 12) of the original series 
are represented by some analytical function F{z;ti,t2)- 
We will further demonstrate how to compute deriva- 
tives F^"^\N) without explicit knowledge of the function 
F{z). We also show how a proper choice of N and n in 

Pj^*^'* allows to control terms neglected in the asymptotic 
expansion. 

In the remainder of this subsection we will implement 
the summation procedure in the form of an efficient nu- 
merical scheme for computing the second-order response 
S'^^^ {ti ,12) at arbitrary times. We reiterate that all terms 
in the series for 5(2) are known in their analytical form 
(via recurrence relations) apart from few parameters de- 
termined by the geometry of the surface, as discussed at 
the end of the previous subsection and in Appendix [Dl 

In the simplest case Nf — 1 the single parameter is 
given by Eq. ([5^ . The series for S'^'^^ is absolutely con- 
vergent, which becomes remarkable only at n > e*^ or 
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^ e*^, and therefore a straightforward computation 
involves an exponentiaUy growing with time number of 
terms with relatively complex structure. At given times 
ti and t2 the terms F{n;ti,t2) of the series constitute 
a sequence of numbers. Our numerical scheme uses the 
procedure described above and reduces the problem to 
the summation of a relatively small and independent of 
time number of terms. This allows to compute S'--^\ti, t2) 
with minimal numerical effort. 

We can justify the numerical procedure separately for 
two overlapping intervals of ti and t2- In the case 
min(ti,i2) ^ 1 the series can be truncated already at 

~ 10. At n > 10 the decay of the terms F{n;ti,t2) is 
exponential, and the remainder is negligible. Therefore, 
the numerical value of S"*-^^ is obtained by the summation 
of several terms. 

In the case min(ii,i2) ^ 1 we are not far from the 
asymptotic region ti,t2 3> 1, and S"*^^^ may be found 
by a simple numerical implementation of the remainder 
calculation ()59|) . Since we use few terms in the expan- 
sion (j59p and approximate them by finite differences, the 
numerical precision is determined by how smoothly the 
numbers F{n; ti, t2) behave as n increases. This is essen- 
tially determined by the value of n. The smoothness is 
only slightly influenced by the values of ii, i2 > 1, which 
can be rationalized by considering the principal asymp- 
totic of 5^^^ at ti,t2 ^ 1. The series ([55]) is not purely 
alternating because of the presence of quasi-oscillations 
(X n^* = e^*l*l'"" in F{n;ti,t2). To ensure a sufhciently 
smooth behavior of F(n; ti, 12) atn>N one must choose 
larger A^ for larger values of \s\. 

In practice, a reasonable relative error < 10""^ is 
achieved if merely A^ ~ 10^ or less terms are retained 
in the partial sum Sn for s ~ 5i. In the expansion ([55]) 
the third term identically vanishes because i?2(l) ~ 0. 
We use two first terms in the expansion, and repre- 
sent the first derivative by the discrete difference with 
the third-order accuracy. Therefore, the neglected con- 
tributions scale at best as 0(fbo(A^; ^i, ^2)6"^™^*!'*^^). 
The efficiency of the procedure is remarkable. We need 
to calculate only around A^ ~ 10^ terms, which should 
be compared with a much bigger number e^*' > 10^ at 
which the summand F{n;ti,t2) starts to decay expo- 
nentially if min(ii,t2) ~ 10. At the first sight, sum- 
mation of 10^ terms versus 10^ does not seem to be a 
crucial improvement, given today's computational capac- 
ities. However, the bottleneck here is computing the in- 
dividual terms that are represented via the hypergeomet- 
ric functions. The computational effort grows dramati- 
cally for the terms F{n;ti,t2) where the series starts to 
converge in the absolute sense. Besides, straightforward 
summation of alternating series is not free from accuracy 
issues. Overall, generating a 2D spectroscopic signal, 
with the developed approach, requires several minutes 
of CPU time using such a simple and high-level tool as 
Mathematica. This qualifies for an "almost analytical" 
calculation of the complete signal. 

In conclusion, we emphasize that only a finite number 



of the series terms is used in the summation procedure. 
Thus, the asymptotic expansion of the second-order re- 
sponse function turns out to be determined by the ex- 
pansions of An{t; s) at t —^ 00. The latter may be viewed 
as spectral decompositions of the evolution operator ma- 
trix elements. In fact, the strongly chaotic system we 
consider is characterized by Ruelle-PoUicott resonances 
i^ii = v + 1/2 ± s. Only even resonances with ly = 2k, 
where k = 0, 1, . . ., contribute to the spectral decompo- 
sition. This result for the second-order response is non- 
trivial because the expansion can be made only after the 
summation of the series over angular harmonics. This 
means that the expansion is only asymptotic. 

The convergence issues can be completely avoided by 
introducing infinitesimal noise. Nonzero noise provides 
the convergence of the series over angular harmonics for 
a given pair of resonances [2^. In the limit of the van- 
ishing diffusion coefficient, the terms of the series turn 
into their noiseless counterparts introduced in Eq. (j56p . 
The application of the remainder summation procedure, 
described above, requires only smoothness of the series 
terms. Therefore, the sum of the series is determined by a 
finite number of terms, and the asymptotic expansion ap- 
pears independent of the vanishing diffusion coefficient. 
This noise regularization is meaningful, since the way the 
convergence is enforced is irrelevant for the asymptotic 
expansion. 

VI. NUMERICAL RESULTS 

Experimental data on 2D time-domain spectroscopy 
that probes the response function S'(^'(ti,t2) is usually 
presented using the so-called 2D spectrum which is ob- 
tained by a numerical 2D Fourier transform of the re- 
sponse function with respect to <i and t2'. 

00 

= JJ dtidt2e'('"i*^+'^^*^)5(2)(ii,t2). (65) 



To provide a spectroscopic view of chaotic vibrational dy- 
namics, in this section we present some numerical results 
on 2D spectra for a particular strongly chaotic system 
studied in this paper. 

As we have shown above, the 2D response in our model 
can be expressed in terms of the properties of the Lapla- 
cian eigenmodes that participate in the expansion of the 
dipole function /(r). Since the dipole is generally a 
smooth function of the system coordinates, the num- 
ber Nf of relevant eigenmodes is typically small. To 
study the important features of the signals we consider 
the second-order response function in the simplest cases 
Nf — 1 and Nf = 2. The latter represents the sim- 
plest situation that involves diagonal and cross peaks in 
the 2D spectrum. Analytical results for this case which 
generalize Section |V] are derived in Appendix iDl 

In the case Nf = 1 the second-order response function 
depends on a spectral parameter s and an additional pa- 
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rameter ag. Both parameters are expressed in terms of 
the only Laplacian eigenmode that represents /(r) [see 
Eqs. jSl), dSl and dUell ]. 

For Nf = 2 we have two spectral parameters si and S2, 
and four additional parameters. Similarly to the Nf = 1 
case, all parameters can be expressed in terms of the 
two relevant Laplacian eigenmodes. As explained at the 
end of Section IV C[ the parameters can be considered as 
independent attributes of a particular Riemann surface. 
For the demonstration of qualitative features, we use a 
particular choice of the parameters. 

The absolute value of the 2D spectroscopic signal 
\ S'''^\lu I, Lu 2)1 is presented in Fig. [3l Panel (a) corre- 
sponds to the Nf = 1 case with the only spectral param- 
eter s — 5i. Panel (b) shows the Nf — 2 case with two 
spectral parameters si = 5i and S2 = 3i. The real and 
imaginary components of the 2D spectra are presented 
in Fig. m In the first case we see a diagonal peak with 
a pronounced stretched feature along coi direction. In 
the second case we also see cross peaks accompanied by 
similar stretched features. 

Diagonal and off-diagonal peaks are also observed in 
spectroscopic signals for harmonic and almost harmonic 
vibrational dynamics. In this case the positions of the 
peaks are given by the frequencies of the underlying pe- 
riodic motions, whereas the width is determined by the 
system-bath interactions and has about the same value 
in both frequency directions. In the chaotic case the peak 
positions are determined by Ruellc-PoUicott resonances, 
rather than by the frequencies of some specific periodic 
orbits. The stretching in spectroscopic signals originates 
from time-domain damped "breathing" oscillations with 
a variable period caused by strong nonlinearity of the un- 
derlying vibrational dynamics and can be interpreted as 
a signature of chaos. 
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FIG. 3: Absolute value of 2D Fourier transform of the second- 
order response function: (a) single resonance s — 5i, (b) linear 
combination of terms with two resonances si = 5i and S2 — 
3i. Linear plots show cross-sections of the spectra at uji — 
0J2 = 5. 







(d) 
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FIG. 4: Real and imaginary parts of 2D spectra: (a) real 
part, (b) imaginary part with the single resonance s = 5i; (c) 
real part, (d) imaginary part of linear combination of terms 
with two different resonances Si = 5i and S2 = 3i. 



VII. CONCLUSION 

In the present paper we studied classical response in a 
strongly chaotic system. Using the Liouville space rep- 
resentation for classical dynamics we found that the re- 
sponse functions exhibit damped oscillations with time. 
The decay is attributed to the mixing property of the 
strongly chaotic dynamics which leads to the efficient 
equilibration in phase space. 

Previous studies of the integrable systems have re- 
vealed the appearance of unphysical divergences in the 
classical nonlinear response functions at long times. The 
divergence originates from the linear time growth of sta- 
bility matrix elements in integrable systems. Chaotic dy- 
namics is characterized by the exponential growth of cer- 
tain stability matrix elements. 

We have described a general qualitative picture of re- 
sponse in strongly mixing (hyperbolic) systems. The ex- 
ponential growth of the stability matrices that can po- 
tentially lead to exponential divergence in response func- 
tions has been interpreted using the Liouville space rep- 
resentation of classical mechanics: As a result of evo- 
lution a smooth initial distribution becomes extremely 
sharp along the stable directions. The second interaction 
with the driving field involves derivatives of the evolved 



17 



distribution along the stable directions. This could re- 
sult in exponentially growing terms in nonlinear response 
functions. We have demonstrated that due to smooth- 
ness of the initial distribution the dangerous terms com- 
pletely cancel out, and the nonlinear response functions 
exhibit exponential decay at long times. Stated differ- 
ently, in mixing systems stability matrices have exponen- 
tially growing and exponentially decaying components. 
Due to the smooth character of the dipole function that 
describes the coupling to the driving field, the growing 
components of the stability matrices are simply elimi- 
nated from the game, and the physical behavior of the 
nonlinear response functions is determined by the expo- 
nentially decaying components. 

To confirm the established qualitative picture we per- 
formed detailed calculation of the linear and second-order 
response for a chaotic model of a free particle moving 
along a compact surface with constant negative curva- 
ture. The model possesses dynamical symmetry that al- 
lows for an exact solution by applying the group rep- 
resentation theory. We found the long-time asymptotic 
behavior of the linear and second-order response func- 
tion that has a form of exponentially damped oscillations. 
Complete asymptotic series obtained in this paper can be 
viewed as expansions in Ruelle-PoUicott resonances. The 
expansion of the linear response in RP resonances is in 
agreement with the earlier results for the two-time cor- 
relation functions [l^l , since the latter is directly related 
to the linear response via the FD theorem. 

The analogue of the RP expansion for the nonlinear re- 
sponse is of more nontrivial nature. The eigenmodes that 
represent the RP resonances are generalized, rather than 
smooth functions. In the linear case the initial smooth 
distribution should be decomposed in the RP modes. The 
signal is computed by convoluting the RP modes with 
the smooth dipole function. Both operations are well- 
defined for generalized functions. In the case of nonlin- 
ear response the second interaction with the driving field 
involves acting on a generalized function with a differ- 
ential operator followed by projecting it onto a general- 
ized function. The legitimacy of the latter operation is 
less obvious, and is related to aforementioned cancella- 
tions of the dangerous terms. The RP decomposition for 
the linear response is represented by a converging series, 
whereas the nonlinear response is given by an asymptotic 
series. Computation of the expansion coefhcients in the 
nonlinear case requires a delicate summation of almost 
sign alternating series. These are other signatures of the 
nontrivial character of the RP decomposition in the non- 
linear case. 

In our forthcoming publication |28| we consider 
Langevin dynamics associated with the model consid- 
ered here. Spectral decomposition in the corresponding 
eigenmodes of the Fokker-Planck operator is a stable le- 
gitimate procedure. The spectral decomposition in the 
presence of noise is a converging series. We show that 
in the limit of vanishing noise the converging series be- 
comes asymptotic and reproduces the expansion derived 



in this paper. The asymptotic expansions of the linear 
and nonlinear response functions can be interpreted as 
decompositions in RP resonances. 

We suggest to apply our results for the interpretation 
of spectroscopic data. Classical chaos is quite generic 
for large molecules. Even for smaller number of nuclear 
degrees of freedom, e.g. in small systems of hydrogen 
bonds [1^, the shape of the effective potential energy 
can make the dynamics chaotic. We have argued that 
restricted motion in the regions with negative curvature 
of the molecular potential can be qualitatively described 
as free motion along a compact Riemann manifold with 
negative curvature. Moreover, the chaos may be caused 
even by the positive but inhomogeneous curvature [30| of 
the configuration space. 

We have considered in detail the second-order response 
which is absent in the spectroscopy of the bulk materials 
or in isotropic environments. Spectroscopy on the sur- 
face [sij . in the absence of the inversion symmetry, can 
measure nonvanishing response functions of even order. 

A dynamical system of the general type has the mixed 
phase space which includes both chaotic regions and sta- 
bility islands. Even in such systems, when the true long- 
time asymptotic decay of correlations is rather power- 
like than exponential, RP resonances may by noticeable, 
leading to the intermediate asymptotic, exponential de- 
cay persisting for a very long time [12] . 

When 2D spectroscopic data is interpreted in terms of 
the underlying dynamics the peaks are usually attributed 
to some periodic motions in the system. We have shown 
that chaotic dynamics can result in similar diagonal and 
off-diagonal peaks in spectroscopic signals that should be 
attributed to RP resonances, and are nor related to any 
specific periodic orbits. Note that the frequencies and 
decay rates of the RP resonances can be retrieved from 
the dynamical ^-function that can be represented as a 
product over all periodic orbits [l3l. Stated differently, 
the RP resonances are related to periodic orbits, yet in an 
extremely collective way, and may not be attributed to 
any particular periodic motions. Therefore, they can be 
referred to as collective chaotic resonances. Our numeri- 
cal results show pronounced stretched features associated 
with diagonal and off-diagonal peaks. These features re- 
sult from breathing damped oscillations that originate 
from contributions of multiple RP resonances with sim- 
ilar oscillation frequencies (imaginary parts) and differ- 
ent damping rates (real parts). Breathing oscillations 
are known to be typical for strongly nonlinear dynamical 
systems. On the other hand, 2D spectra in harmonic or 
almost harmonic systems coupled to a multi-mode har- 
monic bath usually show similar peak patterns in both 
frequency directions. We suggest that the stretched ver- 
sus non-stretched peak shape could have a capacity of 
distinguishing between the collective versus individual 
nature of resonances in spectroscopic signals. 
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APPENDIX A: DYNAMICAL SYMMETRY OF 
GEODESIC FLOWS IN RIEMANN SURFACES 
WITH CONSTANT NEGATIVE CURVATURE 
AND REPRESENTATION THEORY 



This Appendix contains some basic aspects on the ge- 
ometry of Riemann surfaces with constant negative cur- 
vature that are necessary for employing the SO{2, 1) dy- 
namical symmetry. This allows to decompose the free- 
particle dynamics using irreducible representations and 
map the original dynamical problem onto some effective 
1-dimensional dynamics in a circle. 

The fundamental group Tg — Tri{Mg) of a Riemann 
surface of genus g which describes noncontractible 
closed paths is generated by 2*7 elements aj,bj, where 
i — \, . . . ^g, with the only relation 0?=! ~ ^■ 



In the case 5 > 1 the compact surface is covered H Mg 
by the hyperbolic plane that can be implemented as a 
pseudosphere determined by the equation -I- 2/2 ~ 2^0 — 
— 1 with yo > embedded in the 3D Minkowski space 
with a metric di'^ = —dy^ + dy\ -\- (other well-known 
equivalent implementations include the Poincare disc and 
the upper complex half-plane). The fundamental group 
acts freely in H, so that = Tg\H. The reduced phase 
space Mg, considered as a bundle — *■ — M^, 
whose fibers are unit velocity or momentum vectors, be- 
ing pulled back to H forms a bundle ^ G ^ H , where 
G = SO{2, 1) can be represented as a pseudo-orthogonal 
group. This can be visualized as follows: G = 5*0(2, 1) 
consists of points in H (positions) that can be considered 
as 3D vectors with norm —1 in the Minkowski metric, to- 
gether with unit velocity vectors that can be interpreted 
as norm 1 vectors orthogonal to the position vectors (with 
respect to the Minkowski metric). Extending these two 
pseudo-orthonormal vectors to a pseudo-orthonormal ba- 
sis set, we interpret G as the space of pseudo-orthonormal 
basis sets, the latter can be thought of as the pseudo- 
orthogonal group S'0(2, 1). Factorizing G ^ 5*0(2,1) 
with respect to the right action of its maximal compact 
subgroup K ^ 50(2) we arrive at if = G/K. The ac- 
tion of Vg in H can be naturally extended to a left action 
of T g in G, which determines the embedding T g C G. 
This leads to convenient for our purposes representations 
Af2 5^ Tg\G/K and = TgXG. In particular this in- 
terprets the action of G in Mg as originating from right 
action of G in itself. 

The presented picture has a very transparent interpre- 
tation. We start with a Riemann surface Mg of genus 
g > 1 with constant negative scalar curvature K = —1. 
As stated in section IIVI we have three canonical vector 
fields in M^, i.e. ctz = d/d9, cti that determines the 
geodesic flow (classical dynamics of a free particle), and 
o'2 — [ci , CTx] . In the case of constant curvature they form 
the Lie algebra so(2, 1) with respect to the vector field 
commutator. Explicit expressions for the vector fields 
are presented in Section IIVI and in Appendix [B) This 
defines an action of so(2, 1) in M^ that can be natu- 



rally extended to G —i- H, considered as the pull-back of 
Mg Mg to H. The algebra action in G can be inte- 
grated to a group action, which implies that locally G has 
a structure of the universal Lie group associated with the 
algebra so(2, 1). Careful studying of the global properties 
of G shows that in fact G = 50(2, 1). This immediately 
implies that H ^ M^ is equivalent to 50(2, l)/50(2) 
as a Riemann space, i.e. H is the hyperbolic space. This 
has a very important implication that the local structure 
of any Riemann surface Mg with constant scalar curva- 
ture K — —1 is locally equivalent to the hyperbolic space 
H, whereas Mg is locally equivalent to G. In particular, 
this implies that all local quantities such as the Lapla- 
cian operator in Mg, the Poisson bracket ui and the 
geodesic flow in Mg can be actually computed in H and 
G, followed by expressing them using the canonical vec- 
tor fields ai, I = 1,2, z. This is how, e.g. Eq. (|2ip can be 
derived in a straightforward way. 

Finally we note that for a Riemann surface with g > 1 
and not necessarily constant curvature the metric deter- 
mines a complex-analytical structure, and the universal 
cover H — > Alg is equivalent to the hyperbolic plane and 
preserves the complex-analytical structure. This implies 
Afg = rg\H. Since the group of conformal diffeomor- 
phisms of H coincides with its isometry group 50(2, 1), 
we have Tg C G. This defines a metric in Mg that 
has constant curvature K = —1 and is conformally- 
equivalent to the original one. Since that, denoting by 
Mg the moduli space of complex analytical structures for 
genus g we can describe a Riemann surface with constant 
curvature by a pair [rj, K) with r/ G Aig and K < Q being 
the scalar curvature. 

The representations of M^ and A/^ in terms of the 
groups K,T a G 



M^ 



= r\G, 

G/K 



M^ 



T\G/K = T\H = M^/K , (Al) 



are very convenient for describing invariant integration 
measures. Since G, K, and F are unimodular groups, 
according to the theorem on invariant measures in ho- 
mogeneous spaces (23I [23, '25*1 , there is a unique invariant 
measure in F\G that is locally equivalent (in the sense of 
covering) to the Haar invariant measure in G. The latter 
generates an invariant measure in G / K , which is locally 
equivalent to the measure in Af^ generated by the metric. 
Combined with the main property of invariant measures 
in homogeneous spaces, this implies for any integrable 
function g{x) 



dxg{x) 



dr i^g{r,e). (A2) 



This measure provides an invariant scalar product in the 
space H of functions in M^. Besides, implementing func- 
tions g{r) in M'^ as functions g{x) in M^ that do not 
depend on 9, i.e. a^g — 0, we have 



M2 



drg{r) 



M3 



dxgix) , 



(A3) 
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and, therefore Ti can be decomposed into a direct sum of 
irreducible representations of G 

n^n^"'^® H.®0m„7^("). (A4) 

seSpcc(M2) nez 



In Eq. (gH denotes the one-dimensional unit repre- 
sentation that represents constant functions, Tis denote 
the principal series representations (with purely imag- 
inary values of s S Speco(M^) where Ims > 0) and 
complimentary series representations (with real values of 
< s < 1/2), whereas Ti^") are the discrete series rep- 
resentations with the integer factors to„ > describing 
how many times a representation participates in the de- 
composition. 



APPENDIX B: ALGEBRA so(2, 1) AND POISSON 
BRACKET 



In this Appendix we derive the Poisson bracket in the 
form of Eq. (j2ip with the commutation relations between 
the operators ui {I = 1,2, z) given by Eq. ([20]) . 

We start with the canonical form of the Poisson bracket 
in terms of coordinates and conjugated momenta pi = 
dL/dr': 

d d d d 

opi or^ or^ opi 



(Bl) 



where we imply summation over repeated indices. 

Since free motion conserves the absolute value of the 
momentum C,, it is convenient to make another local 
choice of variables in the tangent bundle TM^ . We rep- 
resent the particle momentum in terms of Q and 6, the 
polar angle which determines the momentum direction. 
The algebra element Uz which generates momentum ro- 
tations in TM^ reads = d/dO. 

The third vector field in the reduced phase space is 
described by the differential operator 

+ A (E,.^g^ipi) i^A- - ^ A 

dq''^ ^'\dpkdpn dpndpk 

Next, we show that cti, cr2: '^z form the algebra so(2, 1). 
The coefficients of the first order linear differential op- 
erator [a2,o'z] include the metric tensor and its first 
derivatives. A straightforward yet tedious calculation 
performed in the basis set of four vectors related to 
the canonical variables Pi,q^ results in the relation 
[o'2,o'z] = — CTi. A similar calculation yields the com- 
mutator [(71,172] — —K{r)f7z which contains configura- 
tion space curvature K(r) and is indeed expressed only 
via the remaining operator <Tz- In the case of the con- 
stant negative curvature, rescaling cti and CT2 leads to the 
so(2, 1) commutation relations: 



A number of transformations allow to express the Pois- 
son bracket (|Bip in terms of the operators 9^, cti, (72, Cz. 
The simple form 

d 9 1, , 

W = — (g) 0-1 - CTl (g) — + 7 (0-2 «) - (7^ «) (72) (B3) 

refiects the presence of the dynamical symmetry. 

The first two terms of the Poisson bracket determine 
the phase space velocity. The constant coefficient in front 
of the second term can be alternatively deduced from the 
condition S'(2)(ii,0) = 0. 



APPENDIX C: UNITARY IRREDUCIBLE 
REPRESENTATIONS OF 50(2, 1) 

In this Appendix we describe a convenient implemen- 
tation of unitary irreducible representation of 5*0(2, 1) in 
terms of functions in a circle S^. This provides an afore- 
mentioned mapping of the free-particle dynamics onto a 
1-dimensional problem. Since the group G = 50(2, 1) 
has a double covering SL2[R) 50(2,1), irreducible 
representations of G are provided by even irreducible 
representations of 5^2 (i?), the latter being well-known 
(see e.g., [Ill)- This situation is conceptually close the 
the case of irreducible representations of 50(3) given by 
even, i.e. integer-spin, representations of a double cover 
SU{2) 50(3) of 50(3). We follow the approach of [H 
for 5^2 (i?) and translate it to the language convenient 
for G^s; 50(2,1). 

The space TLs of an irreducible representation of G that 
belongs to the principal or complimentary series has a 
convenient basis set 4'^ of angular harmonics, so that 
(Jz^k = ik'^k- The anti-Hermitian conjugated raising 
and lowering operators a± = <Ji±ia2, Using the allow for 
the the following relations between the basis functions: 



±A: + - - sj ■^k±i , cr,*fe = ik-^k , (CI) 
(*fc,*fcO = 0, forfc^fc'. 



In the case of imaginary s (principal series) ^m('u) 
are normalized, and correspond to normalized functions 
i>k{x;s), i.e. 



Af3 



dx {tpk'{x;s'))* tpk{x;s) = Sk'kSs's 



(C2) 



[cti, 0-2] = CTz , [<yi,az] — (12 , [(J2,<7z 



where dx is the invariant measure in M'^ defined up to a 
constant [23, 24, 25]. We can alternatively view Eq. (|C1[) 
as a definition of a so(2, 1) representation parametrized 
by s. In Appendix [Bl we directly verify that the operators 
a±, Uz satisfy the so(2, 1) commutation relations given by 
Eq. (|20p . The definition of Eq. (|20p becomes clear when 
we implement the representation space Ti^ as a Hilbert 
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vector space of functions ^'(u) in a circle: 



d d 1 - 2s 

tr^ = — , <Ti — sm u- — \ — cos u . 

du du 2 ' 



(C3) 



d 1 - 2s 



(72 = ~ COSU— — 

du 



■ sm u , 



, , . , , . d 1 - 2s 
a±=exp(±zz.)(T*^ + ^- 



^'fc(u) = exp(ifcu) . 



The validity of Eqs. (pO|l and (jCip for the generators cr 
and basis vectors ^fe defined by Eqs. (|C3|) can be verified 
directly. The reason for such a simple representation of 
the generators in the form of first-order differential oper- 
ators is that each representations Tig is induced j23] from 
a one- dimensional representation, parametrized by s, of 
a two-dimensional subgroup AN C G [1^. It is imple- 
mented in the space of functions in the maximal compact 
subgroup K C G where K = 5*0(2) = S^. In the case 
of imaginary s (principal series) the inducing representa- 
tion is unitary. Therefore the induced representation has 
a natural scalar product [23j 



du 
2^ 



{^'{u))*^{u) 



(C4) 



It can be easily verified that the generators ct; for 
Z = 1,2,2; in the implementation of Eq. (|C3[) are anti- 
Hermitian operators with respect to the scalar product 
defined by Eq. 

It also follows from Eq. (|C4|) that the set of vectors '^k 
constitute an orthonormal basis with the positive scalar 
product (cr+^'o, tT+\E'o) = (1 — 4s^)/4 for imaginary s. 
Besides imaginary s, the scalar product is also positively 
defined for real s with — 1 < 2s < 1 (complimentary se- 
ries). In this case the representations Tig are not unitary 
with respect to the scalar product of Eq. (|C4[) . however, 
another scalar product still diagonal in the basis set of 
can be defined (this procedure is also known as rep- 
resentation unitarization [13] )■ We do not use the the 
complimentary series representations in the calculation 
of the response because they lead to the exponential de- 
cay without the oscillatory behavior. 

Irreducible representations 7i*^") of the discrete series 
can be considered as sub-representations Ti^"^ C Hn+i, 
and H^"^ C Hn-i for n > and n < 0, respectively, gen- 
erated by the vectors ^'fc with k > n and k < n, respec- 
tively, which creates a certain inconvenience. This does 
not constitute a major problem, since the discrete repre- 
sentations can be alternatively holomorphically in- 
duced from unitary representations of the maximal com- 
pact subgroup K = 5*0(2). We are not discussing this 
construction here, since discrete representations do not 
contain the zero-momentum state ^Pq, and, therefore, 
they do not contribute to the linear and second-order 
response. 

It is also important to note that irreducible represen- 
tations of the principal and complimentary series charac- 
terized by opposite values of s are unitary equivalent, i.e. 



H-s — Ti-s, which follows from the fact that the represen- 
tation unitarization is determined by the value of s^ [2^ . 
For the principle series representations this property is 
clearly seen from Eq. (|C1[) . The latter implies that in 
the case of imaginary s the vectors {<7±)^ \l/o for opposite 
values of s are different by just phase factors and can be 
connected by a unitary transformation that is diagonal in 
the ^'fc basis set. The purpose of the unitary transforma- 
tion is to compensate the aforementioned phase factors. 
In particular, this justifies the agreement that the spec- 
trum Spec(M^) contains only imaginary s with Ims > 
and real s with < s < 1/2. 

We conclude this appendix by relating the spectrum 
Spec(M^) to the spectrum of the Laplacian in M^. To 
that end we introduce the Casimir operator G that com- 
mutes with all so(2, 1) generators, and, therefore, is a 
constant in any irreducible representation due to the Shur 
lemma 



1 2 

O = --(0-+(T_ +Cr_(T+) + (ct^) , 



(C5) 



[O, ai] = for Z 



l,2,z, 



1 - 4s^ 



for * e 7^, , 



where we used ladder operators a±. We can further in- 
terpret functions in as functions f{x) in inde- 
pendent of the momentum direction 9, i.e. ct^/ = 0. It is 
one of the signatures of the dynamical symmetry that, if 
acting on functions in M^, the Laplacian operator (which 
is proportional to the Hamiltonian of a quantum free par- 
ticle) is expressed in terms of the Casimir operator as 

-y^ = -\{<J+a^ + <J^a+) = C - {a,f , (C6) 
1 _ 4^2 

V tpoix; s) = — i:o{x; s) . 

It follows from Eq. (jC6[) that the spectrum SpeCo(Af^) of 
a Riemann surface is totally determined by the spectrum 
of its Laplacian V^. The spectrum of the Laplacian on a 
compact surface is discrete, and the eigenvalues are neg- 
ative. The functions ipk{x;s) with s € SpeCo(M^) that 
constitute a basis set in the space of relevant distribu- 
tions according to Eqs. ([25| and ([27]) can be expressed 
in terms of the Laplacian eigenfunctions iPq{x] s) using 
ladder operators according to relations (|C1|) . 

For the sake of completeness we note that the functions 
ipk{x; s) for s G SpeCg(M^) and fixed value of fc, known as 
modular forms of degree k (see, e.g. [l^]) are eigenmodes 
of the Laplacian operator still defined by Eq. (|C6p . 
The Laplacian operator can be viewed as the quantum 
Hamiltonian of a free particle moving in homogeneous 
magnetic field, whose intensity is proportional to k. 

Given the function ^'(u) in the circle, one can find its 
counterpart i^{x; s) in representation space Tis- 

oo 2- 

^{x;s)^ j due-'''---^{u)^:u{x-s), (C7) 

A;— — OO Q 

where t/jki^'-, s) are basis functions in TCs- 
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APPENDIX D: MATRIX ELEMENTS IN 
SECOND-ORDER RESPONSE FUNCTION 

This Appendix includes some details necessary 
for the calculation of the second-order response. 
The latter contains integrals of triple products like 
■>p^{x;si){ai^2ipo{x:,S2))i^m{x;s2)- The Appendix de- 
scribes this geometrical part in the calculation of the 
second-order response function, whereas the dynamical 
part is mapped onto the problem in the circle and solved 
Section HVl The anti-Hermitian operators ai {I — 1,2, z) 
play a fundamental role in the description of the dynam- 
ics. They are associated with vector fields that commute 
according to Eqs. (pO|) . 

In Section FyBl we have derived the recurrence relations 
for matrix elements in the simplest case Si = S2 — = s. 
Since ipo is a real- valued function (see Section |TT|, we 
notice that oi = a_i = ao/2. The recurrence relation 
()52|) is symmetric with respect to the sign reversal of k 
and we conclude that a-k = a-k- Therefore, all terms in 
sets afc and are determined by a single real number oq. 
Another way to come to the conclusion is to notice that 



V'*(a;,s) 



r(n + i 



i + .)r(i-s 



-V'-„(a;,s). (Dl) 



Next we calculate the coefficients involved in the 
second-order response in a general situation when the 
coupling / includes contributions from several irreducible 
representations characterized by imaginary numbers si, 
S2, ■ ■ ■ ■ In this case the coefficients are labeled by addi- 
tional indices: 



qrs _ 



dxipl{x;q)'ipo{x;r)->pkix;s) , (D2) 



bl''" = 2 / dx^l{x;q)(T+^po{x;r)i;k-i{x;s) . (D3) 



Choosing the zero momentum eigenfunctions iPq{x; q) to 
be real, ag*^^ remains invariant under permutations in 
{q,r,s}. Employing Eqs. (jCip the recurrence relations 
(|52| can be generalized: 



(2fc + 1 + 2q){2k + 1 - 2s)al''_^^ = 
- (2fc - 1 - 2q){2k - 1 + 2s)a^^^i 
+ (8fc2 + 1 - + - 4s^)ar , 



(D4) 



peaks in 2D spectroscopic signals. The second-order re- 
sponse has a form similar to Eq. ([55)1 : 



(D7) 



1 ^ ^pqr _ („ „Pqr 



n+-+r a-' - \n + - + p a';. 



t2-Q^-nj (^A; {t2 ; p) A„+i (ti ; r) 
t2 + n + 1 j (Al+i{t2;p)An{h;r) 



Since q, r, s adopt only two values si and S2, the coef- 
ficients turns out to be symmetric, ajj^'' — oF^^, and 
Eqs. (|D4p at fc = allow to express a^" via a^" . Since 
the zero harmonics il}o{x; s) are real, the coefficients 
are also real and invariant under any permutation of q, r 
and s. Therefore, all coefficients aj^" are divided into four 
groups and expressed via only four real numbers ag^^^^^ 

The properties of the sequences {a^}"^"^} and {a^""^^^^} 
have been described above. The conditions necessary to 
express all numbers in these sets via the terms at fc = 
are given by af^ = ag*'*/2. 

The third group of matrix elements includes a^^^^*^, 
Q,siS2Si g^jjj a^^*^*^ . We express all these matrix elements 
' using the following relations: 



via On 



1-252 

2(1 - 2si)' 



S1S2S1 _ 1 + 4^2 ^'^l„siais 



2(1-45?) "0 



(D8) 
(D9) 



a 



The fourth group of matrix elements includes a^^''^*^, 
1^'^^'^^ and a^"^*^*^ which obey the relations similar to 
those for the third group but with si and S2 interchanged. 

Thus they expressed in terms of a single real number 

S2S2S1 



a. 







We assume that four numbers d, 



SlSlSi S2S2S2 „SiSiS2 



and 



On 



to be independent. This is rationalized in a 



l/j^^^ = (2fc + 1 - 2q)aj^'' - (2fc + 1 - 2s)a[^^ . (D5) detailed discussion at the end of subsection IVBl 



This is done similarly to the particular case q = r ^ s 
described above. The asymptotic form that generalizes 
Eq. (|E1[) consists of two contributions: 



(D6) 



The numerical results of Section IVII are presented for 
the case of two terms in the expansion ([50)1 of the dipole 
moment /. This is the simplest case that involves cross- 



APPENDIX E: ASYMPTOTIC EXPANSION 

In this appendix we present some details involved in a 
derivation of the asymptotic expansion given by Eq. (j63p . 

First of all the derivation is based on the assumption 
that the the terms F{n) in our series can be represented 
by some smooth function F{z). This is true for the coef- 
ficients An that enter Eq. ([55)1 due to their explicit form 
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given by Eq. (|38|). However, the coefRcients a„ that are 
also a part of Eq. ([55)) are given as a discrete set accord- 
ing to Eq. ([52]) . To show that Ofc can be represented by 
a smooth function a(z) we note that since Uk is a solu- 
tion of a recurrence relation, its A: — > oo behavior can be 
readily analyzed, which results in an expansion 



oo 



EO-r, 
T 



k- 



m=0 



m— 



(El) 



which yields an analytical function a{z) in a form of an 
analytical expansion 



a{z) 



oo 

ar. 



oo 



(E2) 



in powers of l/z at z = oo. 

So far we have demonstrated that the terms F{k) in 
our series are represented by an analytical function F{z), 
and therefore an approach for calculating the asymptotic 
expansion introduced in Section IV CI is in principle legit- 
imate. However, this is not enough yet, since the asymp- 



totic expansion coefRcients S., 



(2) 



involve the derivatives 



d^Fmim2{z) at integer z = k values of the argument, as 
prescribed by Eq. (|62p . Apparently, it requires explicit 
knowledge on the analytical function F{z). However, 
what we really need as a good enough approximation for 
the derivatives in Eq. (|62p so that the neglected terms 
vanish in the limit iV — > oo. This can be achieved by 
implementing a very simple scheme. For an analytical 
function F{z) denote by F^\n) an m-th order approx- 
imation for its derivative d^F{z)\z=N at z ~ N. Then 
the derivatives d^Fmimii^) ^q. l|^ can be safely 
replaced by their approximate values Fm}mim2{N), pro- 
vided m is large enough. The approximate derivatives 
can be linearly expressed in terms of the series terms 



F{n) by solving an obvious system of linear equations 



1=0 



In the remainder of this appendix we show how to con- 
trol the neglected terms in a derivation of the asymptotic 
series. We start with noting that the asymptotic expan- 
sion of Eq. ([63)1 should be understood in a way that for 
any pair (Mi, M2) of positive integers we have 



Ml M2 



mi— m2— 



_|_ Q |^g-2(Afiti + M2t2)~j 



-2(miti+m2t2) 



(E4) 



Full control over the neglected terms can be achieved by 
a proper choice of the intermediate N(ti,t2) as a function 
of ii,t2 ^ 00 and n in Eq. (pT|) . Note that both N{ti,t2) 
and n also depend on Mi and M2. We choose 

7V(ti,t2) -exp(7min(ii,t2)) , (E5) 
so that for small enough 7 we have 1 <C <C 
e*^,e*^ and we can expand F{n;ti,t2) for n < 
N in powers of e"^*^ and e~^*^ according to 
Eq. (|56p . If we keep the terms up to the degree 
(Mi,M2), the neglected terms can be estimated as ~ 
N^iMi+M2)+r (_2 niin(ti, t2)) e-^^^^''^+^^'^\ where 
r — S corresponds to a very loose estimate. By choosing 7 
to be small enough, specifically (2(Mi -I- M2) -I- r) 7 < 2, 
the neglected terms are o (e"^(*^i*i+*^2t2)) ^nd can be 
safely omitted. Finally, no matter how small 7 is we can 
always find M to be large enough (and independent of ti 
and t2), so that 7V-^(ti,t2) = 0(6-2(^^1*1+^^2*2))^ and 
the terms omitted in Eq. ([62|) can be also safely neglected. 
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